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Executive  Summary 


The  Network  Analysis  Connectivity  Model  (NCAM)  is  a  software  package  designed  to  predict 
the  viability  of  communication  links  between  moving  nodes.1  It  consists  of  six  modules:  the 
deployment  module,  the  propagation  module,  the  antenna  module,  the  noise  module,  the  link 
budget  module,  and  the  connectivity  confidence  interval  (CCI)  module.  The  deployment  module 
keeps  track  of  the  platform  latitude,  longitude,  and  elevation  as  a  function  of  time.  The 
propagation  module  uses  the  Terrain  Integrated  Rough  Earth  Model2  program  to  account  for  the 
effect  of  terrain  on  signal  propagation.  The  antenna  module  applies  antenna  properties,  such  as 
gain  and  radiation  pattern,  to  signal  propagation.  The  noise  module  accounts  for  the  effect  of 
radio  noise  from  man-made  sources,  such  as  jamming,  transfonners,  vehicles,  transmission  from 
other  radios,  or  heat  from  the  broadcasting  radio,  and  natural  sources,  such  as  the  background 
cosmic  radiation,  the  Sun,  and  the  Earth.  The  link  budget  module  takes  the  data  from  the 
previous  modules  and  calculates  the  wireless  link  signal-to-noise  ratio  for  each  communications 
link.  The  CCI  module  calculates  the  appropriate  standard  deviation  of  the  path  loss  associated 
with  each  link  and  computes  the  probability  of  successful  signal  reception.  NCAM  uses  the  six 
modules  to  do  the  calculations  necessary  to  perform  the  wireless  network  simulation.  NCAM 
requires  all  six  modules  to  perform  its  calculations. 

The  approach  used  is  to  do  all  the  calculations  needed  to  determine  the  CCI  at  certain  key  times 
called  snapshot  times.  The  computations  at  the  snapshot  times  form  a  series  of  samples  through 
time  of  link  viability,  which,  when  viewed  in  sequence,  shows  how  well  the  nodes  are  linked 
through  time,  much  like  frames  in  a  movie  show  samples  of  cinema  action.  To  predict  the 
viability  of  communication  links  between  moving  nodes  through  time,  it  is  necessary  to  track  the 
location  of  all  nodes. 

The  deployment  module  functions  by  requiring  the  user  to  input  an  outline  of  the  platform’s  path 
using  a  set  of  sample  points  called  waypoints.  Each  waypoint  is  specified  with  the  latitude, 
longitude,  elevation  (read  from  a  Digital  Terrain  Elevation  Data  database  for  a  ground 
platform3),  and  arrival  and  departure  time  of  each  waypoint.  The  deployment  module  then 
interpolates  the  position  of  the  platfonn  at  the  snapshot  times. 


Still,  G.  W.;  Nealon,  J.  F.  The  Case  for  Using  the  Spherical  Model  to  Calculate  the  Interpolated  Points  in  the  Connectivity 
Software  Deployment  Module ;  ARL-TR-4373;  U.S.  Army  Research  Laboratory:  Aberdeen  Proving  Ground,  MD,  February 
2008. 

2Eppink,  D.;  Keubler,  W.;  TIREM/SEM  Handbook',  ECAC-HDBK-93-076;  Department  of  Defense  Electromagnetic 
Compatibility  Analysis  Center:  Annapolis,  MD,  March  1994;  pp  1-134. 

3 

Pablo’s  Mission  Planning.  http://www.mission-planning.com/DTED_Part2.htm  (accessed  28  June  2007). 


For  verification  purposes,  this  report  documents  the  derivation  of  the  equations  that  the 
deployment  module  uses  to  interpolate  the  positions  of  the  moving  communications  platforms 
(or  nodes)  and  the  development  of  the  algorithm  for  adding,  removing,  and  modifying  the 
waypoints.  In  developing  the  equations,  it  was  assumed  that  there  were  two  constant  velocities: 
either  zero  or  some  maximum  cruising  velocity  vcru.  vcru  is  taken  from  a  database,  and  its  value 
is  selected  based  on  terrain  roughness  and  the  platform  type.  When  the  velocity  changes,  it 
increases  or  decreases  at  a  constant  acceleration  rate  of  ±ao.  This  results  in  a  time  vs.  velocity 
graph  that  is  trapezoidal  in  shape.  Furthermore,  the  velocity  of  the  platform  as  it  arrives  and 
leaves  the  waypoint  must  be  the  same.  A  platform  is  allowed  to  change  acceleration  from 
+  (or  -)  ao  to  -  (or  +)  ao  as  it  passes  through  a  waypoint  or  remain  at  zero  if  it  is  traveling  at  vcru. 
The  platform  also  has  the  option  of  lingering  at  a  waypoint. 

The  user  is  given  the  choice  of  specifying  platform  speed  or  arrival  time  at  each  waypoint.  The 
motion  equations  derived  for  interpolation  also  provide  a  lower  and  upper  bound  of  the  node’s 
speed  or  arrival  time  at  each  platfonn.  So,  the  user  is  given  the  maximum  choice  for  specifying 
node  motion  but  is  prevented  from  entering  times  and  speeds  that  are  impossible  for  the  platform 
to  attain. 

The  algorithm  developed  provides  for  the  removal  or  addition  of  waypoints,  changing  the  arrival 
time/velocity  at  waypoints,  or  changing  the  location  of  the  waypoints.  When  the  waypoint 
properties  are  modified,  where  possible,  the  waypoint  velocity  remains  the  same  or  is  modified 
to  be  the  minimum  or  maximum  allowed  for  the  platform,  and  the  arrival/departure  times  are 
recalculated  for  the  subsequent  waypoints. 

The  velocity  profile  was  chosen  so  as  to  make  the  platform  motion  realistic.  A  simple  profile 
would  have  been  unrealistic.  Making  the  profile  more  complicated  would  make  the  motion  more 
realistic,  but  the  derivation  of  the  equations  needed  for  interpolation  would  have  been  more 
involved  than  those  presented  here.  A  more  complex  velocity  profile  would  also  have  required 
more  user  input.  The  velocity  profile  (and  the  subsequent  derivation  of  the  equations  of  motion) 
here  provided  the  best  compromise  between  user  control,  equation  and  algorithm  simplicity,  and 
realistic  platform  motion. 

NCAM  is  being  developed  jointly  by  the  Missile  Defense  Branch  at  Aberdeen  Proving  Ground, 
MD,  and  the  Communications  Electronic  Warfare  Branch  at  Fort  Monmouth,  NJ,  of  the  U.S. 
Anny  Research  Laboratory’s  Survivability/Lethality  Analysis  Directorate.  For  verification 
purposes,  the  details  of  the  algorithm  and  formulae  development  have  been  documented  and 
shown  to  be  consistent  with  the  assumptions  of  the  model’s  description  of  platform  motion. 


viii 


1.  Introduction 


A  new  software  package  called  the  Network  Connectivity  Analysis  Model  (NCAM)  ( 1 ,  pp  3-4) 
is  currently  undergoing  development.  It  is  being  developed  jointly  by  the  Missile  Defense 
Branch  at  Aberdeen  Proving  Ground,  MD,  and  the  Communications  Electronic  Warfare  Branch 
at  Fort  Momnouth,  NJ,  of  the  U.S.  Army  Research  Laboratory’s  Survivability/Lethality  Analysis 
Directorate.  NCAM  predicts  the  viability  of  communications  links  between  moving  nodes. 

Six  modules  constitute  NCAM,  as  illustrated  by  the  screen  capture  of  the  graphical  user  interface 
shown  in  figure  1 .  The  deployment  module  calculates  and  tracks  the  movements  of  the  wireless 
network  nodes.  The  propagation  module  calls  on  the  program  TIREM  (the  Terrain  Integrated 
Rough  Earth  Model)  (2)  to  calculate  signal  attenuation  due  to  the  effect  of  atmosphere  and 
terrain.  The  antenna  module  contains  data  concerning  the  antennae  used,  including  gain,  loss, 
and  electromagnetic  radiation  pattern.  The  noise  module  computes  receiver  noise  attributable  to 
internal,  external,  man-made,  natural,  hostile,  and  nonhostile  noise  sources.  The  link  budget 
module  takes  the  data  from  the  previous  modules  and  calculates  the  wireless  link  signal-to-noise 
(S/N)  ratio  for  each  communications  link.  The  connectivity  confidence  interval  (CCI)  module 
calculates  the  appropriate  standard  deviation  of  the  path  loss  associated  with  each  link  and 
computes  the  probability  of  successful  signal  reception.  NCAM  uses  the  six  modules  to  do  the 
calculations  necessary  to  perform  the  wireless  network  simulation. 


Figure  1.  Modules  of  the  connectivity  software  NCAM. 
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1.1  Background 

The  deployment  module  was  the  first  module  developed.  At  certain  times  in  the  simulation 
called  snapshot  times,  the  S/N  ratio  and  the  CCI  is  calculated  for  every  link  between  the  moving 
nodes.  To  do  those  calculations,  the  distance  between  nodes  must  be  derived  from  their 
positions.  To  obtain  those  distances,  the  Earth  could  have  been  modeled  either  as  a  perfect 
sphere  (with  elevations  in  the  radial  direction,  whose  values  are  very  small  compared  to  the 
radius  of  the  Earth)  or  as  an  oblate  spheroid  (also  with  small  variations  in  the  radial  direction  to 
represent  elevation).  Developers  chose  the  spherical  model  of  the  Earth  as  the  basis  model  for 
distance  computation  because  it  allowed  more  rapid  computation,  with  a  small,  but  acceptable, 
sacrifice  in  accuracy  {1 ,  pp  3-4). 

1.2  Purpose 

The  user  defines  the  path  for  each  communications  node  or  platform  by  selecting  key  locations 
along  the  path  called  waypoints.  Successive  waypoints  in  time  are  then  connected  by  great 
circles,  which  are  the  shortest  lines  between  two  points  on  the  surface  of  a  sphere.  The  great 
circle  approach  accounts  for  the  locations  and  the  elevation  of  each  waypoint. 

This  report  traces  the  derivation  of  equations  describing  platform  motion.  The  equations  derived 
will  be  used  by  NCAM’s  deployment  module  to  obtain  platform  position  as  a  function  of  time. 
This  allows  NCAM  to  calculate  and  use  the  distance  information  in  its  calculation  of  link 
viability.  An  algorithm  allowing  for  the  modification  of  the  number  of  waypoints  for  each 
platform  and  the  modification  of  waypoint  location  is  also  included  to  allow  modification  of  the 
platform  paths.  This  presentation  of  the  motion  equations  and  the  algorithm  is  intended  to 
provide  evidence  for  validating  the  NCAM  software. 

1.3  Scope 

The  derivations  are  limited  to  the  deployment  module.  Because  elevations  are  considered,  the 
equations  derived  can  be  applied  to  either  ground  or  airborne  nodes.  Due  to  the  nature  of  the 
trigonometric  functions,  special  limitations  are  also  considered  with  the  derivation  of  the 
equations. 


2.  Methodology 


Platform  motion  is  assumed  to  be  along  a  great  circle  path.  Platform  velocities  vary  between  a 
zero  or  nonzero  cruise  velocity.  Velocities  between  zero  and  the  cruise  velocity  are  attained 
while  accelerating  between  the  cruise  and  zero  velocity  at  either  a  positive  or  negative  rate  of 
constant  acceleration.  Solutions  to  special  problems  occasioned  by  the  limitations  of  the 
trigonometric  functions  are  also  considered,  as  is  the  process  needed  to  add,  remove,  or  modify 
waypoint  locations. 
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3.  Analysis 


The  analysis  consists  of  deriving  the  equations  of  motion  for  the  platforms  using  spherical 
coordinates.  After  reviewing  the  spherical  coordinate  system,  the  equation  for  a  great  circle  is 
derived.  The  equation  is  derived  for  two  arbitrary  points  that  define  the  circle,  and  special 
circumstances  such  as  the  two  defining  points  (or  waypoints)  nearly  antipodal  on  the  equator,  at 
opposite  poles,  antipodal,  on  the  same  or  opposite  meridians,  one  waypoint  at  the  pole,  both 
waypoints  on  the  equator,  or  one  point  at  the  intersection  of  the  equator,  and  the  international 
date  line,  or  prime  meridian.  The  analysis  ends  with  the  derivation  of  the  value  for  the  Earth’s 
average  radius  used  in  the  spherical  model. 

3.1  Spherical  Latitude-Longitude  Coordinate  System  Review 

The  latitude-longitude  system  of  coordinates  used  in  modern  navigation  on  Earth  is  illustrated  in 
figure  2.  (f>  is  the  latitude,  which  varies  from  -90°  to  90°.  The  designation  “north  latitude” 
indicates  a  positive  value,  while  “south  latitude”  indicates  a  negative  value.  9  is  the  longitude, 
varying  from  -180°  to  +180°.  “East”  corresponds  to  positive  latitude,  while  “west”  connotes 
negative  latitudes.  R  is  the  distance  from  the  Earth’s  center  and  is  the  sum  of  the  mean  sea  level 
and  the  elevation.  R  is  always  positive  and  can  have  a  value  between  0  and  infinity. 


Figure  2.  The  latitude-longitude  coordinate  system  superimposed  on 
the  rectilinear  coordinate  system. 
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To  facilitate  angle  and  distance  calculations,  we  superimpose  a  three-dimensional  (3-D), 
rectilinear  coordinate  system  on  the  latitude  and  longitude  system.  The  x  axis  starts  at  the 
Earth’s  center  point  and  passes  through  the  point  where  the  prime  meridian  (0°  longitude)  and 
equator  (0°  latitude)  intersect.  Likewise,  the  y  axis  starts  at  the  center  point  but  goes  through  the 
point  at  the  intersection  of  the  equator  (0°  latitude)  and  90°  east  longitude.  The  z  axis  begins  at 
the  Earth’s  center  as  well  but  passes  through  the  North  Pole  (90°  north  latitude). 

Converting  from  the  latitude-longitude  coordinate  system  to  the  rectilinear  system  involves 
employing  simple  geometric  relationships.  In  the  vertical  right  triangle  in  figure  2,  the  side 
opposite  the  latitude  angle  cj)  is  the  z  coordinate.  Expressed  in  terms  of  the  latitude-longitude 
system,  this  becomes 

z  =  R  sin  (j)  .  (1) 

The  side  adjacent  to  the  latitude  angle  cj)  becomes  the  hypotenuse  for  the  two  right  triangles  in  the 
x-y  plane  shown  in  figure  2,  with  a  value  of  R  cos  cj).  The  side  opposite  the  longitude  angle  9  is 
the  y  coordinate,  which,  in  tenns  of  R,  cj),  and  9,  is 

y  =  R  cos  cj)  sin  9  .  (2) 

The  side  adjacent  to  the  latitude  angle  cj)  is  the  x  coordinate,  rendered 

x  =  R  cos  cj)  cos  9  .  (3) 

3.2  Equation  of  a  Great  Circle  on  a  Sphere 

We  now  derive  the  equation  of  a  great  circle  on  a  sphere  in  tenns  of  latitude  angle  cj)  and  the 
longitude  angle  9.  First,  consider  the  situation  in  figure  3:  a  sphere  of  arbitrary  radius  Ro  with 
an  equator  and  an  arbitrary  great  circle.  Note  that  the  hidden  lines  are  not  shown  for  clarity.  A 
great  circle  is  a  circle  of  maximum  radius  drawn  on  a  sphere  so  that  it  too  has  a  radius  of  Ro, 
dividing  the  sphere  into  equal  hemispheres  (3).  Note  that  the  equator  is  also  a  great  circle,  and 
that  it,  the  great  circle,  and  the  sphere  share  the  same  center  point.  Included  is  a  rectilinear 
primed  coordinate  system  with  three  axes:  x’,  y’,  and  z’.  The  relationship  between  the 
nonprimed  coordinate  system  and  the  primed  coordinate  system  will  be  established  later.  The  y’ 
axis  begins  at  the  center  point  and  goes  through  one  of  the  two  points  of  intersection  between  the 
equator  and  the  great  circle.  The  x’  axis  begins  at  the  center  point  and  goes  through  a  point  on 
the  equator  such  that  the  x’  axis  forms  a  99°  angle  with  the  y’  axis.  The  z’  axis  begins  at  the 
center  point  and  goes  through  the  pole,  forming  a  right  angle  with  both  the  x’  and  the  y’  axes. 

The  sense  of  the  three  primed  axis  is  such  that  x’  cross  y’  equals  z’. 
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Next,  we  consider  the  angle  the  great  circle  forms  with  the  equator.  Figure  4  shows  the  great 
circle  and  the  equator  with  the  two  intersecting  planes:  plane  1  contains  the  equator  and  the 
sphere’s  center  point,  while  plane  2  contains  the  great  circle  and  the  sphere’s  center  point.  (The 
sphere  has  been  removed  from  figure  4  for  clarity.)  Planes  1  and  2  form  the  angle  (j)max  with  each 
other.  Plane  2  also  forms  the  angle  (j)max  with  the  x’  axis.  The  x’  axis  is  also  in  plane  1,  while  the 
y’  axis  is  in  the  planes’  line  of  intersection.  The  z’  axis  is  normal  to  plane  1.  For  the  next  part  of 
the  derivation,  we  consider  the  appearance  of  the  planes  as  viewed  in  the  direction  of  the  y’  axis, 
along  the  line  of  intersection  of  planes  1  and  2. 

Figure  5  shows  the  planes  from  this  viewpoint.  We  wish  to  concentrate  on  the  relationship  of 
planes  1  and  2  and  realize  the  quantitative  relationship  between  them  and  the  z’  and  x’  axes. 

From  the  point  of  view  of  figure  5,  plane  2  is  a  straight  line  in  the  x’-z’  axis.  So,  in  the  x’-z’ 
system,  the  line  must  obey  the  equation  z’=mx’+b  ( 4 ).  This  is  true  for  every  point  in  plane  2, 
which  contains  the  great  circle.  Since  plane  2  intersects  plane  1  at  the  point  x’  =  z’  =  0,  b 
becomes  zero.  Because  m  is  the  slope  of  plane  2,  this  equals  the  tangent  of  the  angle  4>max.  When 
the  substitutions  established  thus  far  have  been  made,  the  equation  for  a  straight  line  becomes 

z’  =  tan  cj)maX  x’  .  (4) 
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Figure  4.  The  angle  between  the  great  circle  and  the  equator  and  the  planes  containing  them, 
been  removed  for  clarity. 


The  sphere  has 


Figure  5.  Viewing  the  sphere  with  planes  1  and  2  of  figure  4 
along  the  y’  axis.  The  sphere  has  been  removed  for 
clarity. 


Limiting  our  consideration  of  the  equation  for  plane  2  to  the  great  circle  portion  of  plane  2,  we 
substitute  equations  1  and  3  for  z’  and  x’  and  use  the  value  of  Ro  for  the  sphere’s  and  circle’s 
radii  in  planes  1  and  2.  Because  we  are  in  the  primed  coordinate  system  and  the  notation  must  be 
kept  consistent,  the  great  circle  equation  becomes 

Ro  sin  <j)’  =  R0  tan  (j)max  cos  cj)’  cos  9’  .  (5) 

Equation  5  is  the  equation  for  a  great  circle  but  only  when  a  great  circle  intersects  the  equator  in 
the  y’  axis.  To  obtain  the  more  general  equation,  we  superimpose  a  second  set  of  axes  on  the 
primed  system:  the  unprimed  x,  y,  and  z  axes  (figure  6).  The  new  z  axis  is  identical  in 
magnitude  and  direction  to  the  z’  axis.  Both  the  z  and  z’  axes  have  the  same  origin  point.  That 
means  that  z  =  z’.  Hence,  equation  1  shows  that  the  latitude  of  the  primed  and  unprimed  systems 
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are  the  same.  In  other  words,  <(>  =  <t>’.  Since  all  the  primed  and  unprimed  axes  have  the  same 
origin  point,  this  puts  the  new  x  and  y  axes  in  the  same  plane  as  the  old  x’  and  y’  axes  (the 
equator’s  plane).  The  x’  and  y’  axes  are  rotated  from  the  x  and  y  axes  by  the  angle  9o.  To  obtain 
the  equation  in  the  unprimed  coordinate  system,  we  will  view  the  situation  in  the  negative  z  (or 
negative  z’,  since  they  are  the  same)  direction  (see  figure  7). 


Figure  7.  View  of  the  geometry  in  figure  6 
in  the  negative  z  direction. 
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Note  that  plane  2  and  the  great  circle  have  been  removed  for  clarity.  The  x  (and  y)  axis  makes 
an  angle  of  0o  with  the  x’  (and  y’)  axis.  An  arbitrary  angle  0’  in  the  primed  coordinate  system 
(or  0  in  the  unprimed  coordinate  system)  has  been  added.  As  is  evident  from  the  figure,  0  in  the 
unprimed  coordinate  system  equals  0’  +  0o.  Solving  for  0’  and  inserting  the  value  into  equation 
5,  taking  advantage  of  the  fact  that  cj)’  =  cj)  and  dividing  equation  5  by  R0  cos  cj)’,  we  realize  the 
equation  for  a  great  circle: 

tan  (j)  =  tan  ())max  cos  (0  -  0O)  .  (6) 

The  point  on  the  great  circle  that  is  farthest  north  of  the  equator  has  a  longitude  of  0’  =  0°.  This 
follows  from  the  fact  that  the  intersection  of  the  equator  and  the  great  circle  happens  on  the  y’ 
axis,  at  0’  =  90°  and  0’  =  -90°.  Since  the  angle  between  the  planes  that  contains  the  great  circle 
and  equator  is  (j)max,  it  follows  that  this  is  the  latitude  of  the  point  on  the  great  circle  at  0’  =  0. 

The  fact  that  this  is  the  most  northern  point  on  the  great  circle  and,  hence,  the  point  with  the 
greatest  latitude  is  the  reason  it  has  been  designated  (j)max  in  equation  6.  As  is  apparent  from 
equation  6,  the  northern  most  point  on  the  great  circle  with  latitude  (j)max  has  a  longitude  of  0o  in 
the  unprimed  coordinate  system. 

3.3  Establishing  0O  and  <t)max  From  Two  Arbitrary  Points 

Suppose  we  are  given  two  arbitrary  waypoints  with  latitude  and  longitude  coordinates  (4>i,0i)  and 
(cj)2,02).  It  is  possible  for  them  to  define  a  great  circle  on  the  sphere.  This  is  because  the  two 
points  and  the  sphere’s  center  point  define  a  plane.  The  same  plane  intersecting  the  sphere 
defines  a  great  circle.  (Figure  6  shows  how  a  great  circle  is  contained  in  a  plane.  This  is  the 
same  great  circle  shown  in  figure  3.) 

It  is  possible  to  derive  the  equation  for  the  great  circle  that  connects  the  two  points.  By  using 
equation  6  with  (<()i,0i)  and  (cj)2,02),  we  can  obtain  the  values  of  c))max  and  0o.  The  first  step  is  to 
plug  in  (4>i,0i)  and  (cj)2,02)  to  equation  6,  realizing  two  instances  of  equation  6  in  equations  7  and 
8: 


and 


tan  c()i  =  tan  c))max  cos  (0i  -  0O) 


tan  4>2  =  tan  ())max  cos  (02  -  0O)  . 


(7) 

(8) 


We  now  have  two  equations  (7  and  8)  with  two  unknowns:  ())max  and  0o-  Dividing  equation  7  by 
equation  8  results  in  an  expression  that  has  one  of  the  two  unknowns,  0max  (and  tan  0max), 
eliminated: 


tan  (j)i  /  tan  (j)2  =  cos  (0i  -  0o)  /  cos  (02  -  0o)  •  (9) 

Next,  we  use  the  well-known  trigonometric  identity  cos  (A  -  B)  =  cos  A  cos  B  +  sin  A  sin  B  (5). 
After  applying  it  to  the  numerator  and  denominator  of  equation  9,  dividing  that  expression  by  cos 
0o,  then  isolating  0o  and  taking  the  inverse  tangent  (arc  tan),  we  find  the  value  of  0o: 
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9o  =  arc  tan  [(tan  (j)i  cos  02  -  tan  cj)2  cos  0i)  /  (tan  (j)2  sin  0i  -  tan  (j)i  sin  02)]  .  (10) 

There  are  two  solutions  for  0o  such  that  -180°  <  0o  <  180°  because  the  function  tan  0  is  periodic. 
That  is,  the  value  of  tan  9  repeats  itself  every  180°  ( 6)  as  illustrated  in  figure  8.  Because  the 
possible  range  of  longitude  goes  from  -180°  to  +180°  (i.e.,  it  spans  360°),  it  is  possible  to  have 
two  full  periods  of  the  function  tan  9  represented  in  longitude.  Figure  8  shows  this.  As  a  result, 
it  is  possible  to  have  two  values  of  9  for  any  single  value  of  tan  9.  The  primary  value  of  0o  as 
expressed  in  equation  10  will  be  in  either  the  first  or  fourth  quadrants.  This  is  generally  the 
default  value  given  for  the  arc  tan  function  for  most  computational  software  packages  (7).  It  will 
therefore  be  necessary  to  add  a  line  of  code  to  make  sure  that  both  possible  values  of  0o  are 
found.  To  find  the  second  value,  we  subtract  180°  if  equation  10  gives  us  a  value  of  0o  that  is  in 
the  first  quadrant  (i.e.,  between  0°  and  90°)  or  add  180°  if  equation  10  produces  a  value  of  0o  that 
is  in  the  fourth  quadrant  (i.e.,  between  0°  and  -90°). 
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Figure  8.  The  function  tan  9  repeats  itself  every  180°. 

9o  (along  with  t)>2  and  02)  will  be  used  to  determine  the  value  of  4>max-  To  find  cjw,  we  again  use 
equation  6.  Substituting  in  the  values  just  derived  for  0O  with  <j)2  and  02,  dividing  by  cos 
(02  -  9o),  and  taking  the  arc  tan,  we  find 

<t>max  =  arc  tan  [tan  tj)2  /  cos  (02  -  0o)]  •  (11) 

Limiting  the  range  of  latitude  from  -90°  to  +90°,  equation  1 1  gives  one  value  of  c))max  for  any  one 
value  of  9o.  Because  our  use  of  equation  10,  and  adding  or  subtracting  180°  as  required, 
produces  two  values  for  90,  it  follows  that  there  are  two  sets  of  values  of  0o  and  c))max  that  satisfy 
equations  10  and  11. 


9 


Cos  9  is  periodic  in  360°  (5).  To  see  how  having  two  values  of  9o  affects  the  values  of  (j)max  in 
equation  1 1,  we  consider  the  unit  circle  in  figure  9.  The  tenn  cos  (92  -  9o)  suggests  that  we 
examine  how  changing  the  value  of  9o  by  189°  (and  hence  changing  the  value  of  the  equation  1 1 
term  02  -  9o  by  189°)  changes  the  value  of  cj)max.  In  figure  9,  we  see  that  the  red  line,  orange 
lines,  and  black  axis  line  form  two  identical  right  triangles.  So,  magnitude  of  the  value  of  cos  (9 
-  189°)  for  any  value  of  9  is  the  same  as  the  cos  9.  However,  the  side  of  the  triangle 
representing  the  cos  (9  -  189°)  is  on  the  opposite  side  of  the  origin  as  cos  9.  Hence,  the  sign  (but 
not  the  magnitude)  of  the  cosine  function  changes  from  negative  to  positive  or  positive  to 
negative  when  the  argument  9  changes  by  189°. 


Figure  9.  A  unit  circle  showing  the  change  of  cos  9 
when  9  changes  by  180°. 


Therefore,  changing  the  value  of  9o  by  189°  in  equation  1 1  changes  the  sign  of  the  tenn  cos  (02  -  9o) 
but  not  the  magnitude.  Changing  the  sign,  but  not  the  magnitude,  of  the  term  cos  (92  -  9o) 
changes  the  sign  but  not  the  magnitude  of  4>max.  So,  the  two  sets  of  values  for  90  and  (j)max  can  be 
expressed  as  (ctw,9o)  and  (-ctw,9o  ±  189°). 

For  any  point  on  a  great  circle  with  the  coordinates  (cj),9),  the  point  exactly  opposite  it  is  189° 
away  in  longitude  9,  equidistant  from  the  equator  in  latitude  cj),  but  in  the  opposite  hemisphere  so 
it  has  the  opposite  sign.  The  two  points  opposite  each  other  on  a  great  circle  can  be  expressed  as 
(cj),9)  and  (-(j>, 9  ±  189°).  The  points  (c()max,0o)  and  (-  c))max,0o  ±  180°)  are  two  points  opposite  each 
other  on  the  same  great  circle.  So  even  though  there  are  two  sets  of  solutions  to  equations  10  and 
1 1  for  t()max  and  9o,  they  both  define  the  same  unique  great  circle  when  applied  in  equation  6. 

This  is  because  of  the  symmetry  of  the  location  the  two  sets  of  points  on  the  great  circle. 
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Figure  10  illustrates  this  symmetry.  Note  that  part  of  plane  1  has  been  removed  so  that  the  point 
(-<t>max,Oo  ±  180°)  is  visible.  It  is  therefore  possible  for  any  unique  great  circle  described  by 
equation  6  to  have  two  sets  of  values  for  90  and  (j)max.  To  ensure  that  (j)max  is  indeed  the  maximum 
latitude  on  the  great  circle,  we  will  choose  the  value  of  9o  that  makes  c))max  positive.  Furthermore, 
note  that  the  points  ((j>max,9o)  and  (-  (|>max,0o  ±  180°),  the  z/z’  axis,  and  the  x’  axis,  occupy  the 
same  plane.  Figure  1 1  shows  the  relationship  between  the  points  ((|>max,0o)  and  (-  cjw, 

0o  ±  180°),  and  the  x’and  z/z’  axes. 


Figure  11.  The  points  (<t>max,0o),  the  origin,  and  (-  <t>max,9o  ±  180°)  share  the  same  plane  with  the  x’  and  the 
z/z’  axes. 
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3.4  Special  Considerations 

Equation  6  is  the  equation  for  a  great  circle  in  terms  of  the  latitude  (j>  and  longitude  9.  Equations 
10  and  1 1  calculate  the  constants  in  equation  6  from  two  arbitrary  points  on  the  circle,  (4>i,0i) 
and  (c(>2,02).  There  are  some  values  of  (cj>i,9i)  and  (cj)2,02)  that  do  not  allow  equations  10  and  1 1  to 
produce  values  for  (f>max  and  9o.  We  will  now  review  those  values. 

3.4.1  (<|)2,02)  on  the  Equator 

This  means  that  (j)2  =  0.  This  reduces  equation  10  to 

9o  =  arc  tan  (  -1  /  tan  02  )  .  (12) 

To  satisfy  equation  12,  02  and  9o  must  differ  by  ±90°.  Note  that  the  great  circle  intersects  the 
equator  at  two  points.  One  is  at  ((j)2  =  0,02),  as  shown  in  figure  12.  The  other  is  at  (<(>2  =  0, 

02  ±  180°).  (The  second  point  is  obscured  by  the  x’-z/z’  plane  in  figure  12.)  The  point  (c))max,Oo) 
has  the  maximum  latitude  on  the  great  circle.  The  two  points  where  the  great  circle  and  equator 
meet  have  latitude  0°.  So,  by  symmetry,  the  point  with  the  maximum  latitude  (c))max,9o)  must  be 
equidistant  from  the  two  points  where  the  great  circle  and  equator  meet.  That  means 
90  =  02  ±  90°  and  that  (c()max,0o)  must  be  in  the  x’-z/z’  plane,  consistent  with  equation  12. 


To  put  it  quantitatively  and  a  little  more  generally, 

|  0o  -  02  |  =  90°,  or  270°  . 


(13) 
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Due  to  the  periodic  nature  of  the  tangent  function  (see  figure  8),  there  are  two  values  of  9o  that 
satisfy  equation  12  in  the  -180°  to  +180°  range.  Note,  too,  that  both  values  of  0o  satisfy  equation 
13.  Normally,  we  would  use  equation  1 1  and  the  two  values  of  0o  to  find  the  two  corresponding 
values  (j)max  and  take  the  set  of  (j>max  and  0o  that  renders  a  positive  (j)max.  Unfortunately,  cos 
(90 -  02)  =  0  because  9o  -  @2  =  ±90°,  or  ±270°  by  equation  13.  This  causes  the  denominator  of 
equation  1 1  to  be  zero.  Since  (4>2,02)  is  on  the  equator,  (jq  =  tan  (j)2  =  0.  Thus,  equation  1 1  reduces 
to  0  divided  by  0,  rendering  the  value  of  4>max  as  undefined. 

Fortunately,  (4>i,0i)  is  not  on  the  equator.  Equation  1 1  could  have  been  derived  to  express  cjw  in 
tenns  of  0 1  by  rendering  it 

<)>max  =  arc  tan  [tan  (jq  /  cos  (9i  -  0o)]  .  (14) 

Because  (4>i,0i)  is  not  on  the  equator,  (jq  is  not  zero,  and  consequently,  0O  -  9i  does  not  equal 
±90°  or  ±270°.  This  fact  permits  us  to  use  equation  14  instead  of  equation  1 1  to  calculate  the 
two  values  of  (j)max  and  choose  the  positive  value  to  make  the  set  (cj)max,0o). 

Had  it  been  (<|q,0i)  on  the  equator  instead  of  (cj)2,02),  equation  12  would  have  rendered  9o  as 

0o  =  arc  tan  (  -  1  /  tan  9i  )  .  (15) 

Since  equation  1 1  makes  no  use  of  cjq  or  0i,  it  could  have  been  used  to  calculate  the  values  of 
cjw.  Hence,  in  the  instance  of  (4>i,0i)  sitting  on  the  equator  and  ((jq,02)  not  sitting  on  the  equator, 
the  calculation  for  0O  and  cj)max  could  have  been  handled  without  special  consideration. 

3.4.2  (<j>2,02)  or  (4>i,0i)  but  not  Both  at  the  Intersection  of  the  Equator  and  the  International 
Date  Line  or  the  Prime  Meridian 

Should  (<[>2,02)  be  at  the  intersection  of  the  equator  and  the  prime  meridian  or  the  equator  and  the 
international  date  line,  it  would  not  be  possible  to  calculate  the  value  of  0o  using  equation  10. 
Equation  10  is  repeated  here  as  equation  16: 

9o  =  arc  tan  [(tan  cjq  cos  02  -  tan  cjq  cos  0i)  /  (tan  cjq  sin  9i  -  tan  cjq  sin  02)]  .  (16) 

The  reason  this  equation  cannot  be  used  is  that  (jq  =  0°  and  02  =  0°  or  ±180°,  resulting  in 
tan  4>2  =  sin  02  =  0,  causing  the  denominator  to  be  zero  and  the  value  of  0O  to  be  undefined. 

Fortunately,  because  the  equator  and  the  great  circle  meet  at  the  points  ((jq  =  O°,02  =  0°)  and 
(0°,±180°),  then  90  must  be  ±90°.  Once  again,  the  symmetry  argument  is  applied  as  shown  in 
figure  13.  The  figure  illustrates  the  instance  where  ((jq  =  O°,02  =  ±180°).  The  other  point  where 
the  equator  and  the  great  circle  meet  (0°,0°)  is  obscured  by  the  x-z  plane.  To  be  the  point  with 
maximum  latitude,  it  must  be  halfway  between  the  points  where  the  equator  and  the  great  circle 
intersect,  which  is  at  0O  =  ±90°.  With  two  values  of  9o  established,  it  is  possible  to  calculate  2 
values  of  (j)max  using  equation  14,  because  (jq  and  0i  are  not  on  the  equator.  It  is  then  a  matter  of 
selecting  the  value  of  0o  which  resulted  in  a  positive  value  of  (j)max. 
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Alternatively,  should  (4>i,©i)  be  at  the  intersection  of  the  equator  and  the  prime  meridian  or  the 
equator  and  the  international  date  line,  equation  16  becomes  undefined  for  the  same  reason  it  did 
when  ((^2,02)  was  at  the  intersection  of  the  equator  and  the  prime  meridian  (or  the  equator  and  the 
international  date  line).  And  for  the  same  reason,  0o  must  once  again  be  ±90°.  But  this  time, 
equation  1 1  is  used  to  detennine  the  two  values  of  cj)max. 

If  0o  tends  toward  ±90°,  then  tan  0o  tends  toward  ±infinity.  This  is  consistent  with  the 
denominator  of  equation  16  tending  toward  zero.  This  gives  us  confidence  that  in  this  instance, 
our  choice  of  0O  =  ±90°  was  the  correct  one. 

3.4.3  Both  (<))2,02)  and  (<j>i,0i)  on  the  Equator 

If  both  (cj)2,02)  and  (4>i,0i)  are  on  the  equator,  then  the  only  great  circle  that  can  contain  both 
points  is  the  equator.  In  that  instance,  (j)max  =  0°,  and  0O  is  undefined.  The  equation  for  the  great 
circle  is  (j)  =  0°,  and  0  can  be  any  value  between  -180°  and  +180°.  Figure  14  illustrates  this  case. 

3.4.4  Both  (<))2,02)  and  (<j>i,0i)  on  the  Same  or  Opposite  Meridians 

Initially,  we  will  assume  that  at  least  one  of  the  two  points  is  not  on  the  equator  or  at  the  pole. 
This  means  that  both  the  numerator  and  denominator  in  equation  16  will  be  finite  so  that  two 
values  of  0O  can  be  established.  However,  the  value  of  cjwx  cannot  be  found  so  readily.  To  see 
why,  we  consult  figure  15. 

In  this  case,  the  mathematical  relationship  between  0i  and  02  is  that  they  are  either  equal  or  differ 
by  180°.  This  causes  the  plane  of  the  great  circle  (plane  2  in  the  figure)  to  be  at  a  right  angle  to 
the  equatorial  plane.  This  means  that  r|)max  =  90°,  causing  tan  (j)max  to  be  undefined  (see  equation 
14).  Note  that  the  sphere  has  been  removed  from  the  figure  for  clarity. 
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Note,  too,  that  the  value  of  9o  returned  by  equation  16  differs  by  ±90°  from  0i  and  ±90°  from  02, 
which,  in  turn,  is  equal  to  or  ±180°  from  9i.  This  is  consistent  with  the  latitude  9o  sitting 
midway  between  the  latitudes  of  the  intersection  points  of  the  equator  and  the  great  circle.  In 
this  case,  the  great  circle  is  a  meridian  line,  so  the  latitude  of  the  intersection  points  of  the 
meridian  and  the  equator  is  equal  to  @i  and  is  ±180°  from  9i.  However,  due  to  the  symmetry  of 
the  great  circle  and  its  orientation  to  the  equator  (see  figure  15),  it  is  impossible  to  tell  which 
value  of  0o  corresponds  to  the  value  of  4>max  =  90°.  Therefore,  the  equation  of  the  great  circle 
with  both  (4>2,02)  and  (cj)i,0i)  on  the  same  or  opposite  meridians  will  be 

0  =  0i,  while  -90°  <  (j)  <  +90°,  and  if  0i  <  0,  9  also  =  0i  +  180°: 

if  0i  >  0,  0  also  =  0i  -  180°,  while  -90°  <  $  <  +90°.  (17) 

In  equation  17,  there  are  only  two  possible  values  of  0,  yet  (j)  may  take  on  any  value.  Note  that 
equation  17  is  equally  valid  if  (0i,cj)i)  is  replaced  by  (02,cj)2). 

3.4.5  (4>2,02)  or  (<j>i,0i)  at  a  Pole 

This  means  that  cj)i  or  (j)2  equals  ±90°.  In  this  case,  the  equation  for  the  great  circle  becomes 
9  =  0j,  while  -90°  <  cj)  <  +90°  and  if  0,  <  0,  9  also  =  0j  +  180°:  if  0j  >  0,  9  also  =  0j  -  180°, 
while  -90°  <  cj)  <  +90°,  where  j  represents  the  number  of  the  point  that  is  not  at  the  pole.  Note 
that  the  equation  is  very  similar  to  equation  17.  Note  too  that  in  this  case,  9i  may  not  be  equal 
to,  or  differ  from,  02  by  180°.  All  that  is  necessary  is  for  4>i  or  cj)2  to  equal  ±90°. 

3.4.6  (<|)2,02)  and  (<|>i,0i)  Are  Antipodal 

Antipodal  points  are  points  whose  coordinates  are  on  exactly  opposite  sides  of  the  sphere.  This 
means  that  9?  =  9i  ±  180°,  and  cj)2  =  — cj)i .  Applying  this  to  the  definition  for  9o  as  listed  in 
equation  16,  we  find 

9o  =  arc  tan  {[tan  <|q  cos  (0i  ±180°)  -  tan  (-(jq)  cos  0i]  / 

[tan  (— c[)i)  sin  9i  -  tan  <|q  sin  (0i±  180°)]}  .  (18) 

For  the  sine  and  cosine  functions,  changing  the  argument  of  the  functions  by  180°  changes  the 
sign  of  the  function  (see  figure  9).  In  other  words  (9), 

sin  (9  ±  180°)  = -sin  9,  and  cos  (0  ±  180°)  = -cos  0  for  all  9  .  (19) 

This  reduced  equation  18  to 

0o  =  arc  tan  {cos  0i[tan  <|q  +  tan  (-<|q)]  /  sin  0i  [tan  (— 4>i)  +  tan  (jq]}  .  (20) 


16 


The  tangent  function  is  a  periodic,  odd  function  (10)  meaning  that  tan  (-(())  =  -tan  cj).  This  causes 
equation  20  to  reduce  to  zero  divided  by  zero,  independent  of  the  latitude  or  longitude.  Hence, 
we  realize  that  when  the  waypoints  are  antipodal,  they  do  not  define  a  unique  great  circle.  To 
see  why,  we  consult  figure  16. 


We  choose  an  arbitrary  point  (c|)i,9i)  in  plane  1,  which  is  the  radius  distance  R  from  the  center 
point  at  the  intersection  of  the  x,  y,  and  z  axes,  which  is  the  sphere’s  center  point.  Next,  we 
choose  the  antipodal  point  (-(|>i,0i  ±  180°),  also  a  distance  R  from  the  center  point.  Because  the 
antipodal  points  are  opposite  each  other  on  the  sphere,  a  straight  line  can  be  drawn  connecting 
the  antipodal  points  and  the  sphere’s  center  point.  This  straight  line,  2R  in  length,  serves  as  the 
diameter  of  the  circle  in  plane  1  in  the  figure.  In  3-D  space,  however,  one  diameter  can  serve  as 
the  diameter  for  an  infinite  number  of  circles,  all  of  which  are  contained  in  a  sphere  of  radius  R, 
as  shown  by  the  circles  in  planes  2  and  3  in  figure  16.  An  infinite  number  of  great  circles  can  be 
rotated  around  their  common  diameter.  Therefore,  two  antipodal  points  on  a  sphere  do  not 
define  a  great  circle.  Note  too  that  antipodal  points  need  not  be  on  the  equator.  Again,  the 
sphere  and  the  equatorial  plane  have  been  omitted  from  the  figure  for  clarity. 
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If  two  waypoints  are  antipodal  on  the  Earth’s  surface,  the  shortest  distance  between  them  will  be 
on  a  path  drawn  through  the  poles.  This  is  because  the  Earth’s  shape  is  very  slightly  flattened  at 
the  poles  due  to  its  rotation.  So  on  our  perfect  sphere  model,  the  equation  for  the  great  circle  for 
antipodal  points  will  be  the  same  as  for  two  points  on  the  same  or  opposite  meridians  discussed 
in  the  previous  section.  The  equation  is  equation  17,  which  is  repeated  here  for  convenience  as 
equation  2 1 : 

0  =  0i,  while  -90°  <  (j)  <  +90°,  and  if  0i  <  0,  0  also  =  0i  +  180°: 

if  0i  >  0,  0  also  =  0i  -  180°,  while  -90°  <  $  <  +90°.  (21) 

Since  the  points  are  antipodal,  it  is  equidistant  to  go  through  either  the  North  or  South  Pole.  For 
consistency,  we  chose  to  go  through  the  North  Pole  for  all  antipodal  waypoints. 

3.4.7  (<|>2,02)  and  (<|>i,0i)  at  Opposite  Poles 

In  this  case,  cjq  and  (jq  are  at  ±90°.  The  infinite  number  of  great  circles  are  meridians.  In  this 
case,  we  chose  the  equation  that  contains  the  prime  meridian  and  the  international  date  line.  So, 

0  =  0°  and  180°  while  -90°  <  <  90°. 

3.4.8  (c)>2,02)  and  (<J>i,0i)  Nearly  Antipodal  on  the  Equator 

In  this  case,  cjq  =  cjq  =  0,  and  |  02  -  0i|  is  almost,  but  not  quite,  180°.  With  the  waypoints  not 
antipodal,  yet  both  on  the  equator,  the  equation  of  the  great  circle  connecting  them  can  be  found 
using  equations  derived  earlier  in  the  instance  when  both  waypoints  were  on  the  equator.  The 
great  circle  connecting  them  would  be  the  equator. 

This  would  be  true  if  the  Earth  were  a  perfect  sphere.  Unfortunately,  it  is  slightly  flattened  at  the 
poles,  and  when  cartographic  precision  is  needed,  the  Earth  is  modeled  as  an  oblate  spheroid 
(11),  which  is  a  special  case  of  an  ellipsoid  (12),  based  on  the  WGS84  (World  Geodetic  System 
for  1984)  model  (13).  In  this  instance,  it  may  well  be  that  the  shortest  path  that  connects  the 
nearly  antipodal  waypoints  passes  closer  to  the  poles  than  to  the  equator.  By  symmetry,  the  path 
could  go  past  either  the  North  or  the  South  Pole.  For  simplicity,  we  decree  that  the  path  passes 
the  North  Pole  (see  figure  17). 
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y  =  y0.  z  =  0 

(0,  6.)  or  (x0,  y0,  0) 


(0,  e2)  or 
(-  Xq,  y0,  0) 


Figure  17.  The  shortest  path  on  the  oblate  spheroid  connecting  two  equatorial  waypoints 
nearly  opposite  each  other  goes  closer  to  the  poles  than  the  equator.  The  plane 
provides  the  cross  section  of  the  oblate  spheroid. 


The  two  waypoints  can  be  connected  with  a  straight  line  that  goes  through  the  Earth.  Next,  we 
superimpose  a  coordinate  system  such  that  the  x  axis  is  parallel  to  the  line  that  connects  the 
nearly  antipodal  waypoints.  The  line  has  the  equation  x  =  any  real  value,  y  =  yo,  where  yo  is  the 
distance  between  the  line  connecting  the  waypoints  and  the  x  axis,  and  z  =  0.  This  puts  the  line 
connecting  the  waypoints  in  the  x-y  plane.  It  is  also  skew  and  perpendicular  to  the  z  axis. 

The  question  that  arises  is  just  how  close  to  the  pole  the  shortest  path  line  passes  and  what  the 
functional  dependence  of  the  distance  is  on  the  value  of  yo.  To  answer  the  question,  we  examine 
the  cross  section  of  the  oblate  spheroid  in  the  y-z  plane,  viewing  along  the  x  axis  in  the  negative 
x  direction.  Figure  18  shows  the  cross  section.  The  cross  section  taken  through  the  poles  is  an 
ellipse  (14)  with  the  following  equation  (15): 

1  =  y2  /  a2  +  z2  /  c2 ,  (22) 
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Oblate  Spheroid  Cross  Section  £  Shortest  Path  Line 
1  =  y2/  a 2  +  Z2/  C2 


Intersection  Point 


Intersecting  Line  Penetration 
Point  (0,  yQ, 


A 

-*  y 


Figure  18.  The  y-z  cross  section  of  the  oblate  spheroid  of  figure  17. 


where  a  is  the  equatorial  radius  and  c  is  the  polar  radius.  The  shortest  path  line  on  the  surface  of 
the  oblate  spheroid  intersects  the  cross-sectional  ellipse  at  a  point  fairly  near  the  pole.  The 
straight  line  connecting  the  waypoints  intersects  the  y  axis  at  the  point  (0,yo,0).  These  two  points 
are  separated  by  a  distance  D.  The  square  of  the  distance  D  is  represented  by  the  equation 

D2  =  (y-yo)2  +  z2,  (23) 

which  can  be  rendered  in  terms  of  a  single  variable  after  solving  the  elliptical  cross-sectional 
equation  (equation  20)  for  z2  and  substituting  into  equation  23 

D2  =  (y-yo)2  +  c2  (1  -  y2  /  a2)  .  (24) 

It  follows  that  the  shortest  path  line  along  the  surface  of  the  oblate  spheroid  should  be  the  line 
closest  to  the  straight  line  intersecting  the  waypoints.  Therefore,  the  point  on  the  ellipse  where 
the  shortest  mean  sea  level  path  intersects  the  elliptical  cross  section  should  be  a  point  on  the 
ellipse  that  is  closest  to  (0,y0,0).  This  point  should  have  a  value  for  y  such  that  the  distance 
squared,  D  ,  is  minimized.  The  value  for  y  is  easily  found  by  taking  the  y  derivative  of  equation 
24  and  setting  it  to  zero.  The  answer  is 

y  =  y0  /  ( 1  -  c2  /  a2 ) .  (25) 

The  greatest  value  y  can  have  is  the  equatorial  radius  a,  so  substituting  that  value  in  equation  25 
and  solving  for  y0  renders  the  maximum  value  possible  for  y0,  yo-max  rendered 

yo-max  =  a(l -c2/a2) .  (26) 

For  two  points  on  the  equator  and  nearly  antipodal,  the  line  connecting  them  on  the  surface  of  an 
oblate  spheroid  does  not  follow  the  equator  if  they  are  within  two  yo  of  being  180°  opposite  each 
other.  If  they  are  closer  than  that,  the  closest  path  lies  on  the  equator.  Plugging  in  the  WGS84 
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values  for  the  Earth,  for  a  =  6,378,137.0  m  and  c  =  6,356,752.314  245  m  (13),  we  find 
yo-max  =  42697.6727075391  m. 

To  facilitate  the  derivation  of  the  value  of  the  minimum  angle  between  two  waypoints  on  the 
equator  for  which  the  shortest  path  between  them  is  not  on  the  equator,  we  define  the  angle  pmax 
as  follows: 


Pmax  =  arc  sin  (y0-max/  a)  .  (27) 

The  minimum  angle  between  two  waypoints  on  the  equator  for  which  the  shortest  path  is  not  on 
the  equator  then  becomes 

Kmm=  180° -2  Pmax.  (28) 

To  see  why,  we  consult  figure  17.  We  suppose  that  we  are  looking  in  the  negative  z  direction  at 
the  special  case  where  the  waypoints  are  located  on  the  line  y  =  yo-max,  z  =  0,  and  that  we  are 
looking  from  the  North  Pole  (cj)  =  90°)  in  the  negative  z  direction  and  view  the  earth  in  cross 
section  at  the  equator.  We  see  the  layout  as  shown  in  figure  19.  (Note  that  “a”  denotes  the 
equatorial  radius.  The  oblate  spheroid  has  been  removed  for  clarity.) 

Equation  27  describes  the  simple  geometric  relationship  between  yo-max,  the  equatorial  radius  a, 
and  the  angle  pmax  found  in  the  figure.  From  the  figure  is  evident  that  the  angle  Kmm  plus  twice 
Pmax  equals  a  straight  angle,  or  180°.  Solving  for  Kmm  produces  equation  28. 


Figure  19.  Cross-sectional  view  of  the  equator  of  figure  15 

looking  in  the  negative  z  direction  when  y0  is  at  its 
maximum  value  and  the  definition  of  the  angles 
Kmin  and  Pmax  added. 
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Applying  the  values  of  the  WGS84  Geodetic  System  model  (13)  to  equations  27  and  28  gives  us 
the  minimum  separation  for  two  waypoints  on  the  equator  that  requires  the  path  connecting  them 
to  go  off  of  the  equator.  We  conclude  that  if  two  waypoints  on  the  equator  are  between  180°  and 
179.232874830387°  apart  in  longitude,  the  shortest  path  connecting  them  is  not  along  the 
equator.  In  describing  angles  separating  the  waypoints  on  the  equator,  keep  in  mind  the  fact  that 
the  longitude  range  is  between  -180°  and  +180°,  and  that  the  functions  used  to  derive  the  angle 
is  periodic.  Taking  these  facts  into  account,  the  relationship  between  9i  and  02,  which  results  in 
the  shortest  path  between  them  laying  off  of  the  equator,  is 

180°  <  180°  -  1 180°  -  1 02  -  0 1 1 1  <  179.232874830387°.  (29) 

This  is  the  quantitative  definition  of  “nearly  antipodal.” 

The  error  introduced  by  using  a  sphere  to  predict  the  shortest  path  between  two  nearly  antipodal 
points  on  the  equator  can  be  reduced  greatly  by  using  a  special  process.  The  process  involves 
using  the  shortest  path  intersection  point  of  figure  18.  By  symmetry,  we  know  that  the  point  is 
halfway  between  the  two  waypoints  on  the  shortest  path  line  and  midway  on  the  straight  line 
connecting  the  two  waypoints.  (This  is  true  for  any  two  nearly  antipodal  points  on  the  equator — 
the  shortest  path  intersection  point  need  not  be  in  the  y-z  plane,  as  in  figure  17.)  The  situation  is 
illustrated  in  figure  20. 


(o,  e2) 


Straight  Line 
Connecting  (0,  0^) 
Nearly  Antipodal 
Way  Points 


Plane  Perpendicular 
to  Line  Connecting 
Nearly  Antipodal 
Way  Points  which 
Bisects  the  Sphere 


Figure  20.  Two  nearly  antipodal  waypoints  on  the  equator  bisected  by  a  plane  that  is 

perpendicular  to  the  straight  line  connecting  the  nearly  antipodal  waypoints. 
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We  added  a  plane  that  is  defined  by  the  origin,  the  North  and  South  Poles,  and  the  shortest  path 
line  intersection  point  (((>S,0S).  The  plane  is  perpendicular  to  the  straight  line  connecting  the 
nearly  antipodal  points,  and  the  plane  bisects  the  sphere. 

Having  established  that  9S  is  midway  between  9i  and  9?,  we  cannot  simply  sum  9i  and  92  and 
divide  by  2.  The  value  of  9i,  92,  and  9S  is  limited  to  the  interval  -189°  to  +189°.  This  fact,  and 
the  periodic  nature  of  the  variables,  requires  a  more  careful  definition  of  9S.  The  complete 
definition  is  outlined  in  equation  39: 

If  |  92  -  9j  |  <  180°,  9S  =  (9,  +92)/2, 

If  |  92  —  9i  |  >  189°,  and  (  9i  +  92 )  /  2  <  9,  9S  =  [  (  9i  +  92 )  /  2  ]  +  189°, 

If  |  92  -  9i  |  >  189°,  and  (  9!  +  92 )  /  2  >  9,  9s=[(91  +  92)/2]  —  189°.  (39) 

We  now  derive  the  value  of  (j)s,  the  latitude  of  the  shortest  path  line  intersection  point.  We  will 
assume  that  the  relationship  between  the  waypoints,  the  straight  line  intersecting  them,  and  the 
axis  is  as  shown  in  figure  21.  If  the  final  mathematical  expression  defining  cj)s  is  independent  of 
Cartesian  variables,  then  the  expression  will  be  true  for  all  longitudes.  This  is  due  to  the  polar 
axial  symmetry  of  an  oblate  spheroid,  and  the  fact  that  any  ellipse  produced  by  a  cross  section 
that  includes  both  poles  has  the  same  major  and  minor  radii. 


Figure  21.  A  more  general  view  of  figure  17:  the  straight  line 

connecting  the  waypoints  is  not  necessarily  parallel  to 
the  x  or  y  axis,  yet  it  is  perpendicular  to  the  z  axis. 
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To  find  the  value  of  cj)s,  we  first  consider  a  more  general  view  of  figure  17  where  the  line 
connecting  the  waypoints  is  not  necessarily  parallel  to  the  y  or  x  axis,  as  in  figure  21.  Yet  it  is 
still  perpendicular  to  the  z  axis  and  contained  in  the  x-y  plane.  The  first  step  in  finding  (j)s  is  to 
find  the  angle  p .  Keeping  in  mind  the  fact  that  the  limits  of  9i  and  02  are  between  -180°  and 
+  180°,  we  see  from  the  figure  that 

If  |  02  -  0i  |  <  180°,  p  =  (180°-  |  02  -  0!  |  )  /  2  , 

If  |  02-0i  |  >  180°,  p  =  (|  02-0i  |-180°  )/2,  (31) 


and 


h  =  a  sin  p  . 


(32) 


The  next  step  involves  a  more  generalized  view  of  the  cross  section  in  figure  18,  as  shown  in 
figure  22.  The  value  of  h  has  already  been  established  in  equation  32.  We  can  then  find  w  with 
equation  25,  where  yo  is  the  quantity  h,  and  y  becomes  w. 

w  =  h  /  ( 1  -  c2  /  a2)  .  (33) 


In  figure  22,  d  is  the  minimum  distance  between  the  straight  line  connecting  the  waypoints  and 
the  shortest  path  line  intersection  point,  and  r  is  the  distance  from  the  center  point  of  the  oblate 
spheroid  to  the  ellipse.  The  quantity  L  is  found  from  solving  the  equation  for  an  ellipse  using  the 
known  quantities  w,  a,  and  c. 

L2  =  c2  (  1  -  w2  /  a2  )  .  (34) 
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Finally,  the  value  of  (j>s  is  given  as 


(j)s  =  arc  tan  (L  /  w)  .  (35) 

The  final  step  to  approximating  the  path  between  two  nearly  antipodal  waypoints  that  are  on  the 
equator  is  to  construct  two  great  circles:  one  going  from  the  point  (0,9i)  to  (cj)s,0s)  and  one  going 
from  point  ((j)s,9s)  to  (9,92).  The  equation  for  the  two  segments  takes  on  the  form  of  equation  6, 
repeated  here  as  equation  36  for  convenience. 

tan  (f>  =  tan  4>max  cos  (9  -  9o)  .  (36) 

Both  segments  will  have  the  form  of  equation  36,  but  the  values  of  cjw  and  9o  will  be  different. 
For  the  segment  going  from  (9,9i)  to  (<j)s,9s),  the  constants  in  equation  36  will  be  called  cj)max-i 
and  9o-i,  while  the  constants  corresponding  to  the  segment  from  (cj)s,9s)  to  (9,9o)  will  be  called 
<jW-2  and  9o-2-  To  find  the  values,  we  take  advantage  of  the  fact  that  (9,9i)  and  (0,02)  are  on  the 
equator.  Because  of  this,  90-i  and  9o-2  will  be  99°  removed  from(9,9i)  and  (9,92).  So, 

If  |  92  —  9i  |  <  189°  and  92  >  9,,  90_i  =  9,  +  99° 

If  |  92  —  9i  |  <  189°  and  92  <  0h  90-i  =  0i  -  99° 

If  |  92  —  9i  |  >  189°  and  92  >  9h  90-i  =  9j  -  99°  and  if  90-i  <  -189°,  90_i  =  90-i  +  369°. 

If  |  9o  —  9i  |  >  189°  and  ©2  <  9i,  ©o-i  =  9i  +  99°  and  if  ©o-i  >  or  189°,  90-i  =  ©o-i  -  369°.  (37) 

Having  found  the  value  of  ©o-i,  we  can  establish  the  value  of  cj)max-i  by  using  equation  1 1 . 

<t>max-i  =  arc  tan  [tan  (j)s  /  cos  (9S  -  90-i)]  .  (38) 

A  similar  process  is  used  to  find  ©o-2- 

If  I  92  -  9i  |  <  189°  and  ©2  >  ©i,  ©o-2  =  92  -  99° 

If  |  92  -  9i  |  <  189°  and  ©2  <  9,,  90.2  =  92  +  99° 

If  |  ©2  -  9i  |  >  189°  and  02  >  0i,  0O-2  =  02  +  90°  and  if  0O-2  >  180°,  0O-2  =  ©o-i  -  360°. 

If  |  02  -  ©i  |  >  180°  and  02  <  0i,  ©o-2  =  ©2  -  90°  and  if  0o-2  <  or  180°,  0O-2  =  ©o-2  +  360°.  (39) 

As  with  cj)max- 1 ,  max-2  is  found  using  equation  11. 

^max-2  =  arc  tan  [tan  (j)s  /  cos  (0S  -  ©0-2)]  .  (40) 

Note  that  (j)max-i  and  c()max-2  must  always  be  kept  in  the  range  0°  to  90°. 
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The  shortest  line  has  now  been  approximated  as  two  distinct  great  circle  segments.  One  has  the 
equation  tan  (j)  =  tan  (j)max-i  cos  (0  -  9o-i),  for  the  segment  with  0  between  90-i  and  0S,  while  the 
other  has  the  equation  tan  (j)  =  tan  (t)max-2  cos  (9  -  @0-2),  for  the  segment  with  0  between  0O-2  and  0S. 

This  may  appear  to  be  a  crude  approximation  for  the  actual  path  connecting  two  nearly  antipodal 
waypoints  that  lie  on  the  equator.  However,  a  previous  investigation  ( 1 ,  pp  17-19,  27) 
discovered  that  navigation  software  provided  by  the  U.S.  Army  Topographic  Engineering  Center 
called  FORWARD  and  INVERSE  ( 16)  either  could  not  calculate  any  path  between  two  nearly 
antipodal  waypoints  on  a  WGS84  ellipsoid,  or  it  predicted  a  path  along  the  equator 
(7,  pp  17-19,  27).  Although  this  module  uses  a  perfect  sphere  as  its  basis,  it  will  be  able  to 
produce  a  path  closer  to  the  truly  shortest  path  between  nearly  antipodal  points. 

3.4.9  Summary  of  Waypoint  Placement  Problems 

Table  1  summarizes  the  difficulty  of  waypoint  placement  discussed  in  the  previous  sections.  The 
difficulty  arises  when  the  waypoints  are  used  to  calculate  the  parameters  (j)max  and  9o,  which 
define  the  great  circle.  Methods  to  circumvent  their  difficulties  (discussed  in  the  previous 
sections)  are  also  summarized  in  the  table. 

Table  1.  Summary  of  waypoint  placement  problems  and  their  solutions. 


Waypoint  Placement 

Problem 

Solution 

Nearly  antipodal  on  equator. 

Path  connecting  them  is 
nonequatorial;  sphere  model  predicts 
equatorial  path. 

Calculate  intermediate  point  with 
equations  30-35;  use  equations 

36M0  to  define  two  great  circle 
connecting  segments. 

Waypoints  at  opposite  poles. 

Infinite  number  of  meridians  all 
equally  valid. 

Chose  to  use  the  prime  meridian. 

Waypoints  antipodal. 

Antipodal  points  do  not  define  great 
circle. 

Use  great  circle  that  passes  through 
North  Pole. 

Waypoints  on  same  or  opposite 
meridian. 

Equation  1 1  is  undefined. 

Use  great  circle  that  contains 
meridians. 

One  waypoint  at  the  pole. 

Equation  10  is  undefined. 

Use  great  circle  that  contains 
meridians. 

Both  waypoints  on  the  equator. 

Equation  10  is  undefined. 

Use  the  equator  as  great  circle. 

One  waypoint  at  intersection  of 
equator  and  prime  meridian  or 
international  date  line. 

Equation  10  is  undefined. 

Assign  80  =  ±90°,  use  equation  1 1 
or  14  to  calculate  (jw,  use  value  of 
0O,  which  gives  (|>max  >  0. 

3.5  Value  of  the  Earth’s  Average  Radius 

In  a  previous  work  (7,  pp  17-19,  27)  it  was  established  that  modeling  the  Earth  as  a  perfect 
sphere  provides  a  basis  for  interwaypoint  distance  calculation.  The  distances  calculated  would 
be  accurate  enough  to  serve  the  needs  of  the  software.  So,  the  question  arises  as  to  how  to 
establish  the  value  of  the  radius  of  the  perfect  sphere. 
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Our  approach  is  to  pick  a  radius  such  that  if  a  latitude  and  longitude  were  randomly  chosen  on  a 
WGS84  oblate  spheroid,  sea  level  would  have  a  50%  chance  of  being  higher  than  the  average 
radius  and  a  50%  chance  of  being  lower  than  the  average  radius.  This  entails  adding  every 
possible  radius  and  dividing  by  the  number  of  radii.  Since  there  are  an  infinite  number  of  radii  in 
an  oblate  spheroid,  we  integrate  all  possible  radii  over  all  space,  then  divide  by  the  integration 
over  all  space  to  define  the  average  radius  Rave  as  where  r(4>,0)  is  the  radius  of  the  latitude  4>  and 
longitude  9,  and  dS  is  the  solid  angle  differential,  also  a  function  of  latitude  and  longitude. 


Rave  = 


■0=  +180°/ 


'e  =  -180° 


■0=  +180°/ 


'e  =  -180° 


(41) 


The  solid  angle  differential  in  spherical  coordinates  is  written  dS  =  cost))  dc()  d9  (17).  (Note  that 
the  reference’s  differential  is  written  as  dS  =  sin  cj)  dc()  d9.  This  is  because  the  reference’s  angle  cj) 
is  measured  down  from  the  North  Pole.  In  our  coordinate  system,  the  reference  angle  ()>  is 
measured  up  from  the  equator.)  To  obtain  the  radius  of  an  oblate  spheroid,  we  consider  the 
oblate  spheroid’s  general  equation  in  Cartesian  coordinates  (18),  which  is  1  =  (x  +  y  )  /  a  +  z  / 
c  ,  where  a  is  the  equatorial  radius  and  c  is  the  polar  radius.  Rewriting  this  in  spherical 
coordinates  using  equations  1,  2,  and  3,  and  using  r  (the  dummy  variable  in  lieu  of  R),  we  have 

1  =  (  r2  cos  2  9  cos  2  cj)  +  r2  sin  2  9  cos  2  cj)  )  /  a2  +  r2  sin  2  9  /  c2  .  (42) 

When  we  apply  the  trigonometric  identity  cos29  +  sin29  =  1,  solve  for  r,  and  change  the 
integration  limits  to  radians  from  degrees,  using  equation  42,  equation  41  becomes 


R  = 

xave 


•6  =  +7T 


'&  =  -7T 


-  +n 


'<)  =  -7T 


=  +7T/2 

cos  <j>  d <j>  6$ 

cos2  <j>  +  sin2  A 


=  -n/2 


=  -tt/2 


(43) 


Evaluating  the  denominator  is  straightforward:  the  value  is  4n.  Evaluating  the  variable  9  in  the 
numerator  is  also  straightforward:  the  value  is  2n.  This  is  because  an  oblate  spheroid  is 
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symmetric  about  the  z  axis,  making  it  independent  of  0.  With  those  variables  evaluated, 
equation  43  then  becomes 


R  =  — — 

rvave  2 


=  +7T/2 

cos  <s>  6<l> 

cos2  <j>  +  sin2  <j> 


<j>  =  -7t/2 


(44) 


To  make  the  integral  wieldy,  a  variable  substitution  is  made  such  that  u  =  sin  (j).  Taking  the 
differential  du  =  cos  4>  dc(),  changing  the  limits,  and  doing  some  algebraic  manipulation  renders 
the  average  radius  Rave  as 


u  =  1 


R  =  — 

H'ave  2 


ac 


V a2  -  °2 ' 


du 


(45) 


-  cz 


+  uJ 


u  —  -  1 


Evaluating  this  integral  (19,  20),  applying  the  limits  of  integration,  and  perfonning  some 
algebraic  manipulation  render  the  average  radius  as 


R 


ave 


( 


a  +~\/a  -  c  ^ 
a  -~\/a  -  c  / 


—  ac  ArcTanh  "Vi  -  c2  /  a2 

V  a2  -  c2 ' 


(46) 


because  the  equatorial  radius  a  is  greater  than  the  polar  radius  c.  After  applying  the  WGS84 
values  for  the  equatorial  radius  a  and  the  polar  radius  c  as  6,378,137.0  m,  and  c  = 
6,356,752.314245  m,  respectively  (13),  the  value  of  Rave  is  6370994.40182752  m,  which  varies 
by  <0.0008%  from  the  value  offered  by  the  Chemical  Rubber  Handbook  (21),  which  was  not 
derived  using  the  WGS84  oblate  spheroid  model. 


4.  Results 


The  formulae  derived  in  spherical  coordinates  from  the  analysis  of  section  3  is  next  applied  to 
obtaining  the  results  of  algorithms  and  formulae  needed  to  describe  platform  motion.  After  a 
presentation  of  terms  used  to  describe  moving  platforms,  the  equations  of  motion  are  derived 
assuming  constant  acceleration  and  constant  deceleration  up  to  and  down  from  the  platform’s 
cruising  velocity.  Next,  an  algorithm  is  presented  for  adding,  removing,  and  modifying 
waypoints.  Finally,  assuming  that  the  snapshot  times  are  known,  the  location  of  the  platform 
between  waypoints  at  the  snapshot  times  is  obtained. 
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4.1  Definition  of  Moving  Platform  Variables 

For  all  moving  platforms,  we  assume  that  the  velocity  vs.  time  profile  takes  a  consistent  structure 
between  two  waypoints.  The  relevant  variables  and  constants  used  to  describe  the  profile  are 
displayed  in  figure  23.  Asa  platform  approaches  the  jth  waypoint  (labeled  “Waypoint  j:”  in  the 
figure),  it  is  traveling  at  a  constant  maximum  cruise  velocity,  vcru.  This  maximum  cruise  velocity 
is  a  function  of  the  kind  of  platform  being  modeled  (high-mobility  multipurpose  wheeled  vehicle 
[HMMWV]  vs.  tank  vs.  unmanned  aerial  vehicle  [UAV]  vs.  dismounted  Soldier,  etc.)  and  the 
kind  of  terrain  the  platform  is  crossing  (paved  road  vs.  unpaved  road  vs.  smooth  terrain  vs.  rocky 
terrain,  etc.)  assuming  that  it  is  a  ground  platform.  The  value  for  the  vcru  is  read  from  a  database 
of  platform  performance  specifications  provided  by  Army  development  centers.  If  the  platform 
is  part  of  a  squad  or  convoy,  vcru  can  be  assigned  the  value  of  the  maximum  cruising  velocity  of 
the  slowest  vehicle  in  the  squad  or  convoy  to  ensure  that  the  platforms  remain  in  formation. 


Departure  Location  xj+1(<|>j+1,ej+1,Ej+1) 
Departure  Velocity,  vj+1 
Departure  Time,  tdJ+1 


Cruise  at  Constant 
Velocity  vcru 


Initial  Departure 
Acceleration  ag 


Arrival  Acceleration  -  ag 

Cruise  at  Constant 
Velocity  vcru 


Departure  Location  Xj(<))j,0j, Ej) 
Departure  Velocity,  Vj 
Departure  Time,  tdl 


Arrival  Location  xj+1(<|>j+1,ej+1,Ej+1) 
Arrival  Velocity  vj+1 
Arrival  Time  t„ 


Way  Point  j: 

Arrival 
Acceleration  -  a^ 

Cruise  at  Constant. 
Velocity  vcru 


Departure  Acceleration  ag 

Arrival  Location  xj(<|)j,ej,Ej) 
'Arrival  Velocity  Vj 
Arrival  Time  t.. 


Figure  23.  Variables  that  describe  moving  platforms  as  they  travel  between  waypoints. 


Note  that  vcru  is  not  the  maximum  velocity  of  which  the  platform  is  capable;  rather,  it  is  a  typical 
speed  the  platform  obtains  when  it  is  in  travel  mode.  For  example,  for  a  dismounted  Soldier  vcru 
might  be  the  brisk  walking  pace  of  3  mph,  even  though  a  Soldier  doing  double  time  can  travel  at 
nearly  twice  that  rate.  We  assign  the  value  that  is  typical  of  a  dismounted  Soldier  in  travel  mode. 
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When  the  platform  approaches  the  jth  waypoint,  it  begins  to  slow  down  with  an  acceleration  of 
-a0.  The  minus  sign  ensures  that  ao  is  assigned  a  positive  value.  This  value,  like  vcru,  is  also  a 
function  of  the  kind  of  platform  as  well  as  the  terrain  the  platform  crosses.  It,  too,  will  be  read 
from  a  database  with  infonnation  supplied  by  Army  development  centers.  And  as  with  the 
variable  vcm,  this  will  not  be  the  maximum  value  of  which  the  platform  is  capable  but  rather  a 
typical  value  for  the  platform. 

When  the  platform  arrives  at  the  jth  waypoint,  it  will  do  so  with  a  velocity  of  Vj  at  the  arrival 
time  taj.  The  waypoint  location  is  specified  in  latitude  (4>j),  longitude  (0j),  and  elevation  (Ej). 

The  elevation  is  the  distance  above  the  average  sea  level  for  a  perfectly  spherical  Earth,  which  is 
the  average  radius  Rave  from  the  Earth’s  center  point.  The  radius  (notated  as  “R”  in  figure  2)  for 
the  jth  waypoint  is  Rave  +  Ej.  The  latitude  and  longitude  are  specified  by  the  user,  while  the 
elevation  Ej  is  the  elevation  of  the  ground  on  which  the  platform  is  traveling  (from  a  database 
with  information  provided  by  the  Digital  Terrain  Elevation  Data  [DTED]  [22])  plus  the  antenna 
height  above  the  ground.  If  the  platform  is  not  a  ground  platfonn,  the  elevation  Ej  is  specified  by 
the  user. 

When  the  platform  departs  the  jth  waypoint  for  the  j+lth  waypoint,  it  will  begin  to  accelerate 
with  the  value  +ao  at  the  departure  time  tdj.  This  is  the  same  value  with  which  it  decelerated 
before  reaching  the  waypoint,  only  in  the  opposite  direction  of  the  deceleration  -ao.  (The  orange 
hash  in  the  figure  that  crosses  the  blue  path  is  the  point  where  the  platfonn  begins  accelerating  or 
decelerating.)  So  the  value  for  ao  must  be  chosen  to  represent  the  platfonn’ s  normal  acceleration 
and  deceleration.  The  departure  velocity  and  the  arrival  velocity  at  the  jth  waypoint  must  be  the 
same,  Vj.  If  the  platform  does  not  stop  at  the  waypoint  but  simply  passes  through  it,  then  it 
leaves  the  waypoint  with  the  same  speed  with  which  it  arrived.  If  it  does  stop,  then  it  must  leave 
the  waypoint  from  a  stop,  so  the  arrival  and  departure  speeds  are  zero.  If  the  platfonn  does  not 
stop  at  the  waypoint,  then  the  arrival  time  taj  and  departure  time  taj  must  also  be  the  same.  If  the 
platform  does  stop,  then  it  is  possible  for  the  platform  to  loiter  at  the  waypoint  for  a  waiting,  or 
loiter,  time  ty  =  taj  -  taj. 

The  platfonn  accelerates  until  it  reaches  the  cruise  velocity  vcru.  It  then  stops  accelerating, 
traveling  at  a  constant  speed  of  vcm.  When  it  approaches  the  j+lth  waypoint,  it  begins  to 
decelerate  at  -a0,  until  it  reaches  the  j+lth  waypoint  with  speed  Vj+i  at  time  tj+l,a.  As  with  the  jth 
waypoint,  if  the  platform  has  a  speed  of  zero  when  it  reaches  the  j+lth  waypoint,  it  may  loiter  for 
time  ty+i  =  tdj+i  -  taj+i  before  it  leaves  the  waypoint  at  time  taj+i  when  it  accelerates  at  a0  to  reach 
vcru  on  its  way  to  the  j+2th  waypoint. 

The  path  the  platfonn  travels  between  two  waypoints  is  the  most  expeditious  path  between  two 
points  on  a  perfect  sphere:  an  arc  of  a  great  circle  (23).  As  the  platfonn  passes  the  waypoint,  it 
changes  direction  instantaneously  so  as  to  follow  the  great  circle  path  to  the  next  waypoint. 
Because  the  direction  changes,  the  platform’s  velocity  changes  instantly.  Yet,  in  an  attempt  to 
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maintain  a  greater  degree  of  realism  in  modeling  platform  motion,  the  platform’s  speed  is  treated 
as  constant  as  it  passes  through  a  waypoint. 

4.1.1  Mathematical  Description  of  Platform  Motion 

Figure  24  shows  graphs  of  acceleration,  speed,  and  distance  vs.  time  of  platform  motion  between 
waypoints.  To  keep  the  motion  as  simple  as  possible,  we  define  the  acceleration  graph  as  a 
simple  step  function:  a0  at  the  departure  time  for  waypoint  j  taj.  At  time  fi,  it  becomes  0,  and  at 
t2  it  becomes  -ao,  until  it  arrives  at  waypoint  j+1  at  time  taj-+i.  At  all  waypoints,  the  acceleration 
isO. 


Figure  24.  Motion  profile  for  moving  platform  between  waypoints. 


Integrating  the  acceleration  graph  yields  the  speed  vs.  time  graph.  The  speed  begins  at  Vj  at 
waypoint  j,  then  increases  at  a  constant  acceleration  rate  of  ao  until  time  fi  when  the  acceleration 
goes  to  zero,  and  the  speed  reaches  vcru.  The  speed  remains  at  vcm  until  time  t2,  when  it  decreases 
at  a  deceleration  rate  of-ao  to  become  Vj+i  at  time  taj+i. 

To  find  the  distance  covered,  we  integrate  the  speed  graph  in  time.  At  time  fy  the  platform 
moves  away  from  the  waypoint  x(()>j,0j,Ej),  such  that  the  distance  increases  with  time  and  the  rate 
of  distance  change  increases  in  time  until  time  fi.  From  time  fi  to  time  t2,  the  distance  increases 
at  a  constant  rate  in  time.  From  time  t2  until  the  platform  reaches  waypoint  j+1  at  time  taj+i,  the 
distance  is  still  increasing  in  time  but  at  a  decreasing  rate. 
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4.1.2  Mathematical  Relationship  Between  the  Motion-Defining  Variables 

To  begin,  we  define  the  acceleration  of  the  platform  in  terms  of  the  time  variable  t  according  to 
the  graph  in  figure  24. 


At  t  =  tdj,  a(t)  =  0  . 

(47) 

At  tdj  <  t  <  ti,  a(t)  =  a0  . 

(48) 

At  t  =  ti,  a(t)  =  a0  . 

(49) 

At  ti  <  t  <  t2,  a(t)  =  0  . 

(50) 

At  t  =  t2,  a(t)  =  -a0  . 

(51) 

At  t2  <  t  <  tgj+i,  a(t)  =  -ao  • 

(52) 

At  t  =  taj+i,  a(t)  =  0  . 

(53) 

To  find  the  speed  dependence  in  time,  we  integrate  equations  48,  50,  and  52  in  time.  The 
constant  of  integration  will  be  determined  by  making  the  speed  at  times  t<y  and  taj+i  equal  to  v, 
and  Vj+i,  respectively.  We  will  also  use  the  fact  that  the  speed  from  figure  24  shows  that  from 
(and  including)  times  fi  to  t2,  the  speed  is  vcru.  We  will  also  assume  that  the  speed  is  continuous 
for  all  time.  From  the  boundary  condition  we  know  that 

at  t  =  tdj  v(t)  =  vj .  (54) 

Integrating  equation  48  in  time,  we  know  that 

at  tdj  <  t  <  ti,  v(t)  =  a0t  +  C,  (55) 

where  C  is  the  constant  of  integration.  Applying  the  boundary  condition  of  equation  54,  we 
evaluate  C,  and  with  some  algebra  we  find  that 

at  tdj  <  t  <  ti,  v(t)  =  a0  ( t  -tdj)  +  vj .  (56) 

Evaluating  the  expression  at  the  first  inflection  point  ti,  we  see  that  the  speed  should  be  vcru  from 
figure  24,  so 

at  t  =  ti,  v(t)  =  ao  ( ti  -  tdj)  +  vj  =  vcru .  (57) 

Between  the  inflection  times  ti  and  t2,  the  speed  remains  constant  at  vcm,  so  that 

at  ti  <t<t2,  v(t)  =  Vcru,  (58) 
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and  it  follows  from  figure  24  that 


at  t  =  t2,  v(t)  =  vcru ,  (59) 

giving  us  another  boundary  condition.  Integrating  equation  47  in  the  time  window  t2  <  t  <  taj+i, 
we  find  that 


at  t2  <  t  <  taj+i,  v(t)  =  -a0  t  +  C  ,  (60) 

where  applying  equation  59  to  evaluate  the  integration  constant  C  and  rearranging  the  variable 
shows  us  that 


at  t2  <  t  <  taj+1,  v(t)  =  vcru  -  a0  ( t  - 12  )  .  (61) 

To  complete  the  speed  profile  in  time  between  the  waypoints,  we  use  the  speed  at  the  j+lth 
waypoint  to  find  that 

at  t  taj+l,  v(t)  VCru  tl()  (  taj+ 1  —  t2  )  Vj+i  .  (62) 

To  obtain  the  distance  D(t)  traveled  by  the  platform  from  waypoint  j  as  a  function  of  time, 
equations  56,  58,  and  61  are  integrated  in  time.  The  result  for  equation  56’s  integration  is 

at  tdj  <  t  <  ti,  D(t)  =  a0  t2  /  2  -  a0  tdj  t  +  vj  t  +  C  .  (63) 

The  value  of  the  integration  constant  C  is  found  by  knowing  that  D(t)  =  0  at  time  t  =  tdj. 

Solving  for  C  and  rearranging  the  variables,  we  find 

at  tdj  <  t  <  fi,  D(t)  =  a0  ( t2  -  tdj2  )/  2  -  (  a0  tdj  -  v, )  ( t  tdj  )  .  (64) 

The  distance  the  platform  obtains  from  the  waypoint  j  at  the  inflection  time  fi  is  found  by 
applying  fi  to  the  variable  t  in  equation  64,  where 

at  t  =  fi,  D(t)  =  a0  (  fi2  -  tdj2  )/  2  -  (  a0  tdj  -  Vj  )  (  fi  -  tdj  )  .  (65) 


Distance  D(t)  must  increase  from  the  value  outlined  in  equation  65.  It  will  increase  with  the 
constant  speed  vcru  as  shown  in  equation  58.  So,  the  integration  of  equation  58  is  added  to 
equation  65  so  that  in  the  next  time  segment,  where  fi  <  t  <  t2,  the  distance  formula  is  given 

at  fi  <  t  <  t2,  D(t)  =  a0  (  fi2  -  tdj2 )/  2  -  (  a0  tdj  -  vj  )  (  fi  -  tdj )  +  vcni  t  +  C  .  (66) 

To  find  the  integration  constant  C,  we  use  the  fact  that  at  time  fi,  equation  66  equals  equation  64. 
The  result  is  that 

at  fi  <  t  <  t2,  D(t)  =  a0  (  fi2  -  tdj2  )/  2  -  (  a0  tdj  -  v, )  (  fi  -  tdj  )  +  vcru  (t  -  fi  )  .  (67) 

The  distance  of  the  platfonn  D(t)  from  waypoint  j  at  the  inflection  time  t2  is 

at  t  =  t2,  D(t)  =  a0  ( fi2  -  tdj2 )/  2  -  (  a0  tdj  -  vj )  (  fi  -  tdj )  +  vcru  (t2  -  fi  ) .  (68) 
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The  distance  will  increase  from  the  point  where  the  platform  is  located  at  time  t2.  The  rate  at 
which  the  distance  D(t)  will  increase  in  time  is  outlined  in  equation  6 1 ,  so  that  adding  the 
integration  of  equation  6 1  shows  that 

at  t2  <  t  <  taj+1,  D(t)  =  a0  ( tf  -  tdj2  )/  2  -  (  a0  tdj  -  Vj )  ( ti  -  tdJ  )  +  vcru  (t2  -  ti  )  + 

vcru  t  -  a0  ( t2  -  2  t2  t )  /  2  +  C  .  (69) 

Evaluation  of  the  integration  constant  C  is  accomplished  by  equating  equation  69  with  equation 
68  for  time  t  =  t2,  so  the  result  is  that 

at  t2  <  t  <  taj+1,  D(t)  =  a0  ( ti2  -  tdj2  )/  2  -  (  a0  tdJ  —  vj )  ( ti  —  tdj  )  +  vcru  (t2  -  ti  )  + 

vCru  ( t  - 12  )  -  a0  ( t2  - 122  )  /  2  +  a0  t2  ( t  - 12  )  .  (70) 

At  time  t  =  taj-+i,  equation  70  gives  the  distance  between  waypoints  j  and  j+1.  For  this  reason,  we 
notate  the  distance  between  them  as  Dj,  evaluate  equation  70  making  t  =  taj+i,  and  gather  some 
tenns  together  to  simplify  the  equation.  After  rearranging  some  of  the  variables,  we  find  that 

at  t  =  taj+1,  Dj  =  a0  ( t,  -  tdj  )2  /  2  +  Vj  ( t!  -  tdJ  )  +  vcru  (taj+,  -  fi  )  -  a0  (taJ+1  - 12  )2  /  2  .  (71) 

Bringing  back  equation  57  and  equation  62,  and  after  applying  some  algebra  to  isolate  and 
number  the  time  variables,  equations  71-73  are  three  independent  equations  showing  the 
relationship  between  nine  variables. 


tl  -  tdj=  (  Vcru  -  Vj  )  /  a0  . 

(72) 

taj+l  —  t2  —  (  VCru  —  Vj+l  )  /  . 

(73) 

By  combining  equations  71-73,  it  is  possible  to  eliminate  the  inflection  times  fi  and  t2,  which 
delineate  the  times  when  the  platfonn  will  change  speed.  Doing  so  produces  an  equation  that 
defines  the  distance  between  the  jth  and  j+lth  waypoint  Dj  in  tenns  of  the  departure  time  of  the 
jth  waypoint  tdj,  the  arrival  time  at  the  j+lth  waypoint  taj+i,  the  platform  speed  at  the  jth 
waypoint  Vj,  the  platfonn  speed  at  the  j+lth  waypoint  Vj+i,  the  platform  cruising  speed  between 
waypoints  vcru,  and  the  platform  acceleration/deceleration  rate  a0.  This  is  equation  74. 

Dj  =  vcru  (taJ+i  -  tdj)  -  (vcru  -  Vj+i)2  /  (2  a0)  -  (vcru  -  Vj)2  /  (2  a0)  .  (74) 

4.1.3  Restrictions  on  Input  Data  -  General 

When  the  user  specifies  the  location  of  the  j  =  1st  waypoint,  the  latitude  and  longitude  will  need 
to  be  input.  The  latitude  can  range  from  -180°  to  +180°,  while  the  longitude  ranges  from  -90°  to 
+90°.  If  the  waypoint  is  for  a  ground  platform,  the  elevation  is  read  from  software  containing 
elevation  data,  like  the  DTED  database  (22).  If  the  platform  is  not  a  ground  platform,  the  user 
inputs  the  elevation  data,  restricted  by  the  ground  elevation,  and  the  platform’s  service  ceiling, 
which  is  read  from  a  database  of  platform  specifications.  By  default,  the  corresponding  arrival 
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time  ta,i  will  be  set  to  0  s,  since  this  is  the  beginning  of  the  simulation.  Also  by  default,  the 
platform  speed  will  be  set  to  0  m/s  if  it  is  a  ground  platform  or  the  minimum  air  speed  if  it  is  an 
air  platform.  The  minimum  air  speed  will  be  read  from  a  database  with  specification  data  for  the 
platform.  The  user  will  then  set  the  loiter  time,  defining  the  departure  time  for  the  first  waypoint 
td,i,  which,  for  j  =  1,  will  be  equal  to  the  loiter  time.  The  cruise  speed  vcra  and  the  acceleration  a0 
are  read  from  the  platform  specification  database. 

Next,  the  user  inputs  the  latitude  and  longitude  for  the  second  waypoint.  The  elevation  is  also 
input  if  the  platform  is  an  air  platfonn,  otherwise  the  elevation  is  read  from  a  DTED-like 
database.  This  defines  the  distance  between  the  waypoints  Di.  Consulting  equation  74  for  j  =  1, 
all  variables  but  two  have  been  defined:  ta,2  and  V2.  If  the  user  defines  the  time  ta,2,  then 
waypoint  speed  V2  is  also  defined.  If  the  user  defines  the  waypoint  speed  V2,  then  by  equation 
74,  ta  2  is  also  defined.  If  the  user  specifies  a  loiter  time  at  waypoint  two  (i.e.,  making  it  so  that 
t3j2  is  not  equal  to  td,2),  then  V2  must  be  equal  to  0  m/s.  The  result  in  that  case  is  that  ta.2  will  be 
calculated  according  to  equation  74,  independent  of  any  user  desire. 

In  the  general  case  (i.e.,  for  j  >1),  the  distance  Dj  is  defined  when  the  user  inputs  the  waypoint 
latitudes,  longitudes,  and  elevations  (or  the  elevations  are  read  from  an  elevation  database).  vcru 
is  also  read  from  a  platfonn  specification  database,  as  is  a0.  tj  and  Vj  were  defined  when  the 
calculations  for  the  last  waypoints  were  done,  leaving  the  user  free  to  choose  either  taj+i,  which 
defines  Vj+i  by  equation  74,  or  Vj+i,  which  defines  taj+i.  Should  the  user  decide  that  the  platform 
loiter  at  the  j+lth  waypoint,  then  Vj+i  must  be  zero,  defining  the  arrival  time  taj+i  and  giving  the 
user  no  flexibility  insofar  as  having  variables  to  define. 

Equation  74  is  rewritten  to  show  the  arrival  time  and  the  speed  at  the  j+lth  waypoint  as  a 
function  of  the  other  variables.  Equation  73  was  obtained  by  simply  solving  equation  75  for 
taj+l- 

taj+i  =  tdj  +  [  2  a0  Dj  +  (vcru  -  Vj+i)2  +  (vcru  -  Vj)2  ]  /  (2  a0  vcru  )  .  (75) 

Equation  76  was  derived  by  solving  equation  74  for  the  speed  at  the  j+lth  waypoint.  When  the 
variable  Vj+i  was  isolated  in  equation  74,  it  became  evident  that  there  were  two  possible  solutions 
for  Vj+i,  as  it  was  necessary  to  take  a  square  root  to  find  Vj+i.  The  two  solutions  found  were  vcru 
plus  an  expression  under  the  radical  and  vcm  minus  an  expression  under  the  radical.  We  take  the 
minus  case  because  the  radical’s  value  is  always  positive,  since  a  platform ’s  speed  cannot  exceed 
Vcru-  Once  the  speed  at  the  j+lth  waypoint  is  calculated  with  equation  76, 

VJ+1  =  vcru  -  o/  2  a0  [  vcm  (taj+1  -  tdJ )  -  Dj  ]  -  ( vcru  -  ^  )2  ,  (76) 

and  the  j+lth  waypoint  arrival  time  taj+i  is  established,  equations  72  and  73  can  be  used  to  find 
the  inflection  times  fi  and  t2. 

When  the  last  waypoint  is  assigned,  we  assume  that  the  speed  of  the  platform  becomes  zero  at 
the  final  waypoint.  Hence,  the  user  has  the  freedom  to  assign  the  location  of  the  final  waypoint, 
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while  the  arrival  time  is  determined  by  equation  75  for  the  final  point.  Should  the  simulation 
stop  before  the  platform  reaches  the  final  waypoint,  the  final  time  is  taken  as  an  interpolation 
point,  and  the  final  latitude,  longitude,  and  elevation  are  calculated,  treating  the  final  time  like  an 
interpolation  time.  Should  the  simulation  end  time  occur  after  the  time  the  platfonn  reaches  the 
final  waypoint,  the  platform  will  simply  treat  the  time  between  the  arrival  time  at  the  final 
waypoint  taj+i  and  the  simulation  end  time  as  a  loiter  time  at  the  waypoint. 

4.1.4  Distance  Calculation 

Using  the  spherical  model,  we  next  calculate  the  distance  between  two  points  notated 
(Rave+Ej,(|)j,0j)  and  (RaVe+Ej+i,cj)j+i,0j+i).  Rave  is  the  average  radius  of  the  Earth;  Ej,  (j)j,  and  0j  are 
the  jth’s  waypoint  elevation,  latitude,  and  longitude;  and  EJ+i,  c)>j+i,  and  0j+i  are  the  j+lth’s 
waypoint  elevation,  latitude,  and  longitude.  The  points  are  shown  with  the  great  circle  they 
define,  and  a  wedge  containing  them  and  the  center  point  in  figure  25.  The  equator  is  in  blue,  the 
great  circle  in  red,  and  the  wedge  fonned  by  the  waypoints  and  the  Earth’s  center  point  is  shown 
in  green.  Note  that  the  waypoints  are  not  on  the  great  circle  but  are  elevated  above  it  by  Ej  and 
Ej+i,  respectively. 


Figure  25.  Two  waypoints  and  the  Earth’s  center  point 
form  a  wedge. 
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To  find  the  distance  between  the  waypoints,  we  must  remove  the  sphere  and  take  a  vantage  point 
normal  to  the  wedge  formed  by  the  two  waypoints  and  the  Earth’s  center.  This  is  shown  in 
figure  26.  Rave  is  the  Earth’s  average  sea  level  radius,  Ej  and  Ej+i  are  the  jth  and  j+lth  waypoint 
elevations,  cj)j  and  (j)j+i  are  the  latitudes  of  the  waypoints,  and  0j  and  9j+i  are  the  longitudes.  A0  is 
the  angle  fonned  at  the  Earth’s  center  point  between  the  lines  connecting  the  waypoints  with  the 
center  point.  Dj  is  the  distance  along  the  platform’s  path  between  the  jth  and  j+lth  waypoint, 
while  A  is  the  line-of-sight  distance  between  the  waypoint  projections  on  the  sea  level  surface  of 
the  Earth. 


First,  we  find  the  relationship  between  the  angle  A9  and  the  line-of-sight  distance  A.  From  the 
two  back-to-back  right  triangles,  we  see  that 

A  /  Rave  =  2  sin  (  A0  /  2  )  .  (77) 

Written  another  way,  equation  77  becomes 

A9  =  2  arc  sin  [  A  /  (  2  Rave  )  ]  .  (78) 

Rave  is  known  but  A  is  not.  To  find  A  we  employ  the  well-known  distance  formula  in  a  3-D, 
rectilinear  coordinate  system  (24). 

A  =  A/  (Xj+1  -  Xj)2  +  (yj+1  -  y,)2  +  (zJ+1  -  Zj)2  .  (79) 

Substituting  equations  1,  2,  and  3,  into  79,  then  employing  equation  78,  the  term  Rave  cancels  to 
reveal  A9  in  terms  of  the  latitudes  and  longitudes  of  the  two  waypoints. 

2 

A9  =  2  arc  sin  {  [(  1  /  2  )  [(cos  cj)j+i  cos  0j+i  -  cos4>j  cos  0j)“  + 

2  2  1/2 

(cos  t()j+i  sin  0j+i  -  cos  ())j  sin  0j)  +  (sin  c()j+i  -  sin  ())j)"]  }  .  (80) 
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Note  that  the  sign  of  A9  is  positive  because  we  are  taking  the  positive  root  of  equation  78. 

Next,  to  find  the  arc  length  Dj  we  employ  the  definition  of  arc  length  in  planer  motion  (25), 
which  is  appropriate  to  figure  26  because  the  wedge  is  contained  in  a  single  plane. 


We  rewrite  equation  81  so  that  the  distance  becomes  S,  a  continuous  variable  as  a  function  of  the 
continuous  variable  9.  We  keep  the  variable  R  to  be  the  distance  from  the  Earth’s  center  point  to 
the  platform.  So,  for  any  point  on  the  path  S  the  platform  traverses,  the  angle  it  makes  with  the 
previous  waypoint  9,  as  measured  at  the  Earth’s  center  point,  is 


We  know  that  the  elevation  at  a  waypoint  Ej  is  related  to  the  distance  from  the  Earth’s  center 
point  Rj  in  that 

Rj  =  Rave  +  Ej  .  (83) 

We  now  use  the  variables  E  and  R  to  represent  the  continuous  elevation  and  continuous  distance 
of  the  platform  from  the  Earth’s  center  point  on  the  path  between  waypoints.  Updating  equation 
83,  we  see  that 

R  =  Rave  +  E  .  (84) 

Differentiating  equation  84  in  time,  it  becomes  self-evident  that  the  rate  of  change  of  elevation 
dE/dt  is  equal  to  the  rate  of  change  of  radius  dR/dt,  as  stated  in  equation  83. 

dR/dt=dE/dt.  (85) 

Now  we  establish  the  relationship  between  the  change  in  time  of  path  length  S  and  the  change  in 
time  of  the  distance  of  the  platfonn  from  the  Earth’s  center  point  R.  For  an  aerial  platform  like  a 
fixed  wing  aircraft  or  a  helicopter,  we  assume  that  the  rate  of  climb  or  descent  is  directly 
proportional  to  its  airspeed  along  the  path  S.  For  a  ground  platform,  the  equivalent  assumption  is 
that  the  elevation  change  is  proportional  to  the  distance  traversed.  Put  another  way,  we  assume 
that  the  slope  of  the  ground  between  waypoints  is  constant.  Expressed  quantitatively,  in  both  the 
aerial  platform  and  the  ground  platform  case, 
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dE  /  dt  =  dR  /  dt  =  k  dS  /  dt , 


(86) 


where  k  is  a  constant  whose  value  will  be  determined  in  due  course. 

Next,  we  wish  to  establish  the  relationship  between  the  platfonn’s  distance  from  the  Earth’s 
center  R,  the  path  length  between  waypoints  S,  and  the  angle  the  platform  makes  at  the  Earth’s 
center  between  the  platform’s  location  and  the  previous  or  jth  waypoint.  Multiplying  equation 
86  by  the  time  differential  dt  and  integrating,  we  see  that 

R  =  k  S  +  C  ,  (87) 

where  C  is  the  constant  of  integration.  To  evaluate  it,  we  know  that  when  S  =  0,  the  platform  is 
at  the  jth  waypoint,  so  the  distance  from  the  Earth’s  center  to  the  waypoint  will  be  Rj  =  Rave  +  Ej. 
This  allows  us  to  evaluate  the  integration  constant  in  equation  87.  After  the  constant’s 
evaluation,  equation  87  is  rendered 

R  =  kS  +  Rj=kS  +  Rave  +  Ej .  (88) 

To  evaluate  k,  it  will  be  necessary  to  differentiate  equation  88  by  9. 

dR  /  d9  =  k  dS  /  d9  .  (89) 

Note  that  the  differential  of  the  constant  terms  is  zero.  Next,  we  define  the  antideferential  F(9) 
such  that  dF(0)  /  d0  =  [  (  dR  /  d9  )2  +  R2  ] 1/2.  Rewriting  equation  82  in  terms  of  F(9),  we  arrive  at 
the  distance  traveled  from  the  last  waypoint  to  the  platform. 

S  =  F(0)  -  F(0)  .  (90) 

Differentiating  equation  90  and  substituting  for  the  term  dS  /  d9  from  equation  89  (keeping  in 
mind  that  F(0)  is  a  constant,  so  its  derivative  is  zero),  we  find 

dR/d0  =  k  [  (  dR  /  d0  )2  +  R2]1/2.  (91) 

Squaring  equation  9 1  and  grouping  the  differentials  together,  then  taking  the  square  root,  we 
arrive  at  a  conventional  first-order  differential  equation. 

(dR  /  R)  =  ±  d0  k  /  (1  -  k2)1/2 .  (92) 

As  the  platform  moves,  the  value  of  9  increases.  If  the  elevation  increases,  then  the  value  of  k 
must  be  positive  in  equation  92.  If  the  elevation  decreases,  then  the  value  of  k  must  be  negative 
in  equation  92,  since  we  assume  that  the  quantity  (1  -  k")  "is  always  positive.  Consider  the 
constant  k  as  listed  in  equation  88.  If  the  value  of  the  elevation  of  the  platform  increases  with 
distance,  k  is  positive.  Likewise,  as  it  decreases  with  distance,  k  is  negative.  So,  we  will  simply 
let  the  sign  of  k  change  to  indicate  the  elevation  change  of  the  platform:  if  k  is  positive,  the 
platform  is  ascending.  If  k  is  negative,  the  platform  is  descending.  To  maintain  consistency  in 
equation  92,  we  take  the  positive  root. 
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(93) 


Solving  equation  92  and  taking  the  natural  antilogarithm,  the  equation  becomes 

R  =  C  exp  [6  k  /  ( 1  -  k2)1/2] , 

where  R  is  the  distance  of  the  platform  to  the  Earth’s  center  point,  9  is  the  angular  distance  along 
the  path,  exp  [  ]  signifies  the  value  e  (roughly  2.718281828459  [26])  raised  to  the  power  of  the 
expression  inside  the  brackets,  and  C  is  the  constant  of  integration.  To  find  C,  we  use  initial 
conditions:  when  9  =  0,  the  right  half  of  equation  93  becomes  C,  which  is  equal  to  the  left  half, 
which  is  R  =  Ri  =  Rave  +  Et. 

To  find  the  value  of  k,  we  evaluate  equation  93  when  the  platfonn  reaches  the  j+ 1th  point  and 
equation  80  to  find  the  angular  distance  between  the  waypoints. 

Rj+1  =  Rave  +  E  j+i=  Rj  exp  [A9  k  /  ( 1  -  k2)1/2]  =  (Rave  +  Ej  )  exp  [AO  k  /  ( 1  -  k2)1/2]  .  (94) 

Taking  the  natural  logarithm  of  equation  94  and  solving  for  k,  we  find  that 

k  =  in  [  (Rave  +  Ej+1)  /  (Rave  +  Ej)  ]  /  {  AO2  +  {  in  [(Rave  +  EJ+1)  /  (Rave  +  Ej)]}2  } 1/2 ,  (95) 

transforming  k  to  a  known  value.  Next,  we  incorporate  equation  88  and  use  the  values  of  the 
variables  when  the  platfonn  reaches  the  j+  1th  waypoint. 

Rj+i  =  Rave  +  Ej+i  =  k  Dj  +  Rj  =  k  Dj  +  Rave  +  Ej .  (96) 

After  solving  for  the  distance  the  platform  travels  between  the  waypoints  Dj  and  eliminating  the 
average  radius  of  the  Earth,  we  find  the  distance  in  terms  of  k  and  the  elevations  at  each 
waypoint. 

Dj  =  (Ej+i  -  Ej)  /  k  .  (97) 

Combining  equations  94  and  96  and  eliminating  Rj+i,  the  radius  at  the  j+lth  waypoint,  and 
algebraically  reananging  some  variables,  it  is  possible  to  find  an  alternative  expression  for  the 
distance  the  platfonn  travels. 

Dj  =  (Rj  /  k )  {exp  [AO  k  /  (1  -  k2)1/2]  -  1}  = 

[  (Rave  +  Ej  )/  k  ]  {exp  [AO  k  /  ( 1  -  k2)172]  -  1 }  .  (98) 

The  constant  k  denotes  the  change  of  elevation  divided  by  the  distance  traveled.  So,  the  possible 
ranges  of  k  can  be  expressed  as  -1  <  k  <  1.  Equation  98  becomes  undefined  if  k  =  ±1,  but 
equation  97  will  still  produce  a  finite  value  for  Dj.  However,  both  equations  97  and  98  are 
undefined  if  k  =  0.  It  is  possible  to  take  the  limit  as  k  approaches  zero  in  equation  98  and  use  the 
result  to  compute  Dj.  We  find 

Dj  =  Rj  AO  =  (Rave  +  Ej)  AO  (99) 
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consistent  with  the  fonnula  for  the  length  of  a  circle’s  arc  (27) — the  radius  times  the  angle  in 
radians.  Since  k  =  0,  the  platform  does  not  change  altitude,  making  the  path  the  arc  of  a  circle. 
Having  established  the  platfonn  path  distance  Dj,  it  is  now  possible  to  calculate  the  arrival  time 
taj+i  using  equation  75  or  the  arrival  velocity  Vj+i  using  equation  76. 

4.1.5  Restrictions  on  Input  Data  -  Maximum  Speed,  Minimum  Time 

The  maximum  speed  a  platfonn  can  attain  between  waypoints  is  vcru,  the  cruising  speed.  To 
arrive  at  the  j+lth  waypoint  in  the  minimum  amount  of  time,  the  speed  profile  would  have  to  be 
as  shown  in  figure  27.  The  profile  differs  from  the  one  displayed  in  figure  24  only  in  that  the 
final  speed  remains  at  the  cruising  velocity,  thus  minimizing  the  arrival  time  taj+i .  Additionally, 
the  two  inflection  times  fi  and  t2  have  been  replaced  with  a  single  inflection  time  f.  It  follows 
that  the  arrival  speed  Vj+i  =  vcru.  Since  figure  27  is  a  special  case  of  figure  24,  the  equation 
describing  the  arrival  time  taj+i  is  a  special  case  of  equation  75.  Equation  75  is  restated  here  as 
equation  100. 

taj+1  =  tdj  +  [  2  a0  Dj  +  (vcru  -  Vj+i)2  +  (Vcru  -  vj)2  ]  /  (2  a0  vcru  )  .  (1 00) 

To  find  the  minimum  arrival  time,  vcru  is  set  equal  to  vJ+i  in  equation  100,  and  the  arrival  time  is 
notated  as  tmin. 

tmin  =  tdj  +  [  2  a0  Dj  +  (vcm  -  Vj)2  ]  /  (2  a0  vcru  )  .  (101) 


For  some  speeds  at  waypoint  j  and  some  distances  between  waypoints  j  and  j+1,  the  platform 
may  not  have  time  to  accelerate  to  the  speed  vcru.  In  that  case,  the  motion  profile  for  speed  vs. 
time  looks  like  figure  28.  The  distance  covered  Dj  is  the  average  speed  (vj  +  vmax)/2  times  the 
time  to  travel  between  waypoints  tmm  -  tdj,  or 

Dj  =  (  tmin  -  tdj  )  (Vj  +  Vmax)  /  2  .  (102) 
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The  relationship  between  the  time  traveled  between  waypoints  tmm  -  tdj,  the  acceleration  a0,  and 
the  speeds  Vj  and  vmax  is 

(  Vmax  —  Vj  )  =  a0  (  tmin  -  tdj  )  .  (103) 

After  eliminating  the  variables  tmm  and  tdj  using  equations  102  and  103,  then  solving  for  vmax, 
we  find  the  maximum  speed  in  this  case  is 

vmax  =  (  2  a0  Dj  +  Vj2)1/2.  (104) 

Using  equation  104  to  eliminate  vmax  from  equation  103,  and  solving  for  tmin,  we  find  the 
corresponding  minimum  arrival  time  at  the  j+ 1th  waypoint  is 

tmin  =  tdj  +  [(2  a0  Dj  +  Vj2) 1/2  -  Vj  ]  /  a0 .  (105) 

So,  the  maximum  speed  at  the  j+ 1th  waypoint  Vj+i  the  user  can  select  for  the  platform  is  vmax,  as 
defined  by  equation  104.  If  the  value  of  equation  104  does  not  exceed  vcru,  then  the 
corresponding  minimum  arrival  time  taj+i  the  user  may  select  is  tmin,  defined  by  equation  105.  If 
the  value  for  vmax  as  defined  by  equation  104  exceeds  the  platform  cruising  speed  vcm,  then  the 
maximum  arrival  speed  the  user  can  select  for  Vj+i  is  vcm  with  a  minimum  arrival  time  taj+i,  for 
the  platform  is  defined  by  equation  105. 

4.1.6  Restrictions  on  Input  Data  -  Minimum  Speed,  Maximum  Time 

The  minimum  speed  with  which  a  platform  can  arrive  at  the  j+ 1th  waypoint  is  generally  (but  not 
always)  zero.  This  is  a  special  case  of  the  velocity  profile  of  figure  24,  shown  here  as  figure  29. 
Equations  75  and  76  govern  the  platform’s  motion  (with  equations  72  and  73  defining  the 
inflection  times  ti  and  t2)  with  Vj+i  set  to  zero.  The  corresponding  maximum  time  is 


tmaX  tdj  T  [  2  ao  Dj  vcru~  "I-  (vcru  Vj)  ]  /  (2  ao  vcru  ) 

(106) 

and  the  inflection  times  are 

1 1  —  tdj  T  (  vcra  —  Vj  )  /  ao 

(107) 

and 

t2  —  tmax  —  ao  /  V  cru  . 

(108) 
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4.1.7  Restrictions  on  Input  Data  -  Peak  Speed  Less  Than  vcru. 

For  some  speeds  at  waypoint  j  and  some  distances  between  waypoints  j  and  j+1,  the  platform 
may  not  be  able  to  achieve  the  cruising  speed  vcru.  In  that  case,  the  speed  profde  resembles 
figure  30.  To  find  the  value  of  vpk,  we  invoke  equations  71-73.  But  in  this  case,  there  is  only 
one  inflection  point  f,  and  at  no  point  between  the  jth  and  j+lth  waypoints  does  the  platform 
travel  at  a  constant  speed.  The  arrival  speed  Vj+i  is  zero.  With  these  considerations  in  mind, 
equation  7 1  for  the  distance  can  be  rewritten  as 

Dj  =  (Vpk  +  Vj)  (ti  -  td j)  /  2  +  vpk  (taJ+i  -  tj)  /  2  .  (109) 


Equations  72  and  73  relate  the  inflection  times  and  the  start  and  ending  times  with  the 
acceleration  and  speed.  They  are 

ti  —  tdj  =  (  Vpk  -  Vj  )  /  a0  (110) 

and 


ta,j+i  —  L — Vpk/-ao.  (Ill) 

We  use  equations  110  and  1 1 1  to  eliminate  the  departure  time  tdj,  the  inflection  time  t,,  and  the 
arrival  time  taj+i  from  equation  109,  rendering  it  as 

Dj  =  (Vpk2  -  Vj2)  /  (2  ao)  +  Vpk2  /  (2  a0)  .  (112) 
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Upon  solving  for  vpk,  we  find 


vpk  =  [  (2  a0  Dj  +  Vj2)  /  2  ] 1/2 .  (113) 

To  find  the  corresponding  maximum  time  tmax,  we  eliminate  the  inflection  time  t,  in  equations 
110  and  111. 


td,j  +  (  Vpk  -  Vj  )  /  a0  =  taj+1  -  Vpk  /  a0 .  (114) 

The  final  step  is  to  substitute  equation  115  for  vpk  into  equation  1 16  and  solve  for  taj+i  and  define 
it  as  tmax. 

tmax  =  tdj+{2[(2a0Dj  +  vj2)/2]1/2-vj}/a0.  (115) 

4.1.8  Restrictions  on  Input  Data  -  Insufficient  Time  for  Platform  to  Reach  Zero 

For  some  speeds  at  waypoint  j  and  some  distances  between  waypoints  j  and  j+1,  the  platform 
may  not  have  time  to  come  to  a  full  stop.  In  that  case,  the  motion  profile  for  speed  vs.  time  looks 
like  figure  31.  In  this  case,  the  arrival  velocity  is  the  minimum  velocity  calculated  with  equation 
116  but  with  the  sign  for  the  acceleration  a0  changed. 

Vmin  =  (Vj2  -  2  a0  Dj)1/2.  (116) 


Likewise,  the  expression  for  the  maximum  arrival  time  resembles  equation  105,  but  the  sign  for 
the  acceleration  is  changed. 

W  =  tdj  +  [Vj  -  (Vj2  -  2  a0  Dj  )1/2]  /  a0 .  (117) 

So,  the  minimum  speed  at  the  j+lth  waypoint  Vj+i  the  user  can  select  for  the  platform  is  zero,  but 
only  if  the  quantity  Vj  -  2  ao  Dj  (see  equation  109)  is  negative.  The  corresponding  maximum 
time  tmax  is  defined  by  equation  106,  with  107  and  108  giving  the  values  for  the  inflection  times 
ti  and  t2.  If  the  quantity  Vj  -  2  ao  Dj  is  positive,  the  value  for  the  minimum  arrival  time  is  set  by 
equation  109.  The  corresponding  maximum  time  tmax  is  calculated  with  equation  117. 


44 


Should  the  user  wish  to  have  a  loiter  time,  the  arrival  speed  va,j+i  must  be  set  to  zero  by  default. 
However,  should  the  user  desire  a  loiter  time  but  the  condition  that  the  quantity  vf  -  2  ao  Dj  is 
>0,  then  the  arrival  speed  is  set  according  to  equation  117,  and  the  loiter  time  is  automatically  set 
to  zero. 

4.1.9  Inputting  Speed  or  Time  Data 

To  facilitate  computing  the  data  needed  for  interpolation,  we  first  consider  the  speed  profile  as 
shown  in  figure  32.  It  defines  the  top  speed  vtop  of  the  platfonn  between  waypoints  for  the 
profile  shown  in  the  figure,  which  happens  at  the  inflection  time  tt.  The  figure  is  a  special  case 
of  the  speed  profile  of  figure  24.  Likewise,  the  quantitative  description  of  the  figure  is  a  special 
case  of  equations  72,  73,  75,  and  76.  To  make  use  of  those  equations,  the  expression  for  vtop 
must  be  derived. 


Suppose  the  user  specifies  the  arrival  velocity  Vj+i.  The  distance  between  the  waypoints  will  be 

Dj  =  (Vtop  +  Vj)  (ti  -  tdj)  /  2  +  (Vj+1  +  Vtop)  (taj+1  —  ti)  /  2  .  (118) 

The  inflection  time  t,  and  start  and  ending  times  tdj  and  taj+i  are  related  to  the  acceleration  and 
velocities  Vj,  Vj+i  by  equations  1 19  and  120. 

tt  -  tdj  =  (  vtop  -Vj)/a0.  (119) 

taj+1  -  tt  =  (  Vj+i  -  Vtop  )  /  -ao .  (120) 

After  substituting  equations  119  and  120  into  equation  1 18  to  eliminate  all  the  time  variables,  we 
find 


Dj  =  (vtop2  -  Vj2)  /  (  2  a0  )  -  (vj+i2  -  Vtop2)  /  (  2  a0  )  . 


Solving  for  vtop  renders 


vtop  -  [  (2  ao  Dj  +  Vj+i2  +  Vj"  )  /  2  ] 1/_  . 


(121) 

(122) 
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If  the  user  specifies  the  arrival  time  instead,  we  first  use  equations  119  and  120  to  eliminate  tt 
and  solve  for  Vj+i. 


Vj+1  2  Vtop  Vj  3-0  (  taj+1  tdj  )  . 

(123) 

We  next  solve  equation  1 15  for  Vj+i. 

v  j+1  =  (2  Vpk2  -  Vj2  -  2  a0  Dj  )1/2 . 

(124) 

Finally,  we  equate  equations  123  and  124  to  eliminate  Vj+i.  After  solving  for  vtop,  we  see  that 

Vtop  =  Vj  +  a0  (taj+1  -  tdj  )  -  [  Vj  a0  (taj+1  -  tdj  )  -  a0  Dj  +  (  1  /  2  )  a02  ( taJ+1  -  tdJ  f  ] 1/2 .  (125) 

When  solving  to  find  equation  125,  the  negative  case  was  selected  on  the  basis  of  calculating  test 
values  for  the  variables. 

Lastly,  we  restate  equations  72,  73,  75,  and  76  with  vtop  substituted  for  vcra,  and  isolate  the 


inflection  times  ti  and  t2. 

ti  =  tdj +  (  vtoP -Vj)/a0.  (126) 

t2  =  taj+1  -  (  Vtop  -  Vj+i  )  /  a0 .  (127) 

taj+i  =  t<ij  +  [  2  a0  Dj  +  (vtop  -  Vj+i)2  +  (vtop  -  vj)2  ]  /  (2  a0  vtop  )  .  (128) 

Vj+i  =  vtop-  {2  a0  [  vt0p  (taj+i  -  tdj)  -  Dj  ]  -  (vtop  -  Vj)2  } 1/2 .  (129) 


If  the  user  selects  the  velocity  Vj+i,  we  calculate  the  top  speed  vtop  using  equation  122.  If  vtop  is 
less  than  vcru,  vtop  retains  its  value.  If  vtop  is  greater  than  vcra,  then  vtop  is  set  equal  to  vcm.  We 
then  use  equation  128  to  calculate  taj+i,  and  equations  126  and  127  to  calculate  the  inflection 
point  times  fi  and  t2. 

If  the  user  decides  to  select  the  j+  1th  arrival  time  taj+i,  we  calculate  the  top  speed  vtop  using 
equation  125.  If  vtop  is  less  than  vcru,  vtop  retains  its  value.  If  vtop  is  greater  than  vcru,  then  vtop  is 
set  equal  to  vcru.  We  then  use  equation  129  to  calculate  Vj+i,  and  equations  126  and  127  to 
calculate  the  inflection  point  times  fi  and  t2. 

4,1,10  Summary  of  Inputting  Waypoint  Data 

Figure  33  shows  the  overall  process  of  inputting  waypoint  data.  The  process  shown  in  the  figure 
is  inclusive  of  equations  47  to  129,  although  it  does  not  include  the  restrictions  outlined  in  table 
1 .  The  rectangular  boxes  represent  processes  to  be  executed.  The  diamond  boxes  are  questions 
with  yes  or  no  answers.  The  answers  detennine  the  calculations,  tasks,  and  the  order  in  which 
the  tasks  and  calculations  are  perfonned.  The  green  letters  indicate  the  “Start”  box,  and  the 
variables  defining  the  platform’s  motion  are  in  blue.  Violet  shows  the  indices,  amber  indicates 
equation  numbers,  and  red  shows  when  it  is  necessary  to  “Do  Interpolation.” 

At  the  beginning,  the  user  identifies  the  platform  type  and  the  number  of  waypoints  (the 
recording  index),  while  the  software  retrieves  the  acceleration  rate  a0  and  the  cruising  speed  vcru 
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Figure  33.  Overview  of  the  waypoint  data  input  process. 
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based  on  the  platform  type.  If  the  platfonn  is  not  a  ground  platform,  the  user  inputs  latitude, 
longitude,  and  elevation.  In  the  case  that  the  platform  is  a  ground  platform,  the  software  reads 
the  elevation  data  from  the  DTED  database  (22).  The  first  point  motion  data,  Vi  and  ta,i,  are  set 
to  zero.  The  user  inputs  the  loiter  time,  and  the  index  j  is  set  to  1.  Implicit  in  these  processes  is 
the  assumption  that  every  platform  has  at  least  two  distinct  waypoints.  The  case  where  a 
platform  has  only  one  waypoint  will  be  treated  later. 

If  the  index  indicates  that  the  next  waypoint  is  the  last,  a  special  process  begins.  As  with  the  first 
waypoint,  the  user  inputs  the  latitude,  longitude,  and  (if  the  platfonn  is  not  a  ground  platform) 
elevation.  When  this  data  is  used,  the  distance  from  the  jth  waypoint  is  calculated.  If  the  distance 
is  not  sufficient  to  allow  the  platfonn  to  stop  (vf  is  not  >2  ao  Dj  [remember,  this  is  the  final 
waypoint  and  the  platform  must  anive  here  with  a  speed  of  zero])  while  it  decelerates  at  ao,  the 
deployment  module  assigns  the  minimum  possible  speed  at  the  waypoint  Vj+i.  It  then  calculates 
the  conesponding  anival  time  taj+i,  departure  time  tci,ri  (remember,  there  is  no  loiter  time  because 
the  speed  is  not  zero),  and  inflection  time  f.  Lastly,  it  forces  the  user  to  input  another  waypoint  (in 
this  case,  fi  =  t2,  so  the  i  represents  the  subscripts  1  and  2  since  there  is  only  one  inflection  point). 
This  process  continues  until  the  platform  traverses  enough  distance  to  come  to  a  halt.  If  there  is 
enough  distance  for  the  platform  to  come  to  a  halt,  then  the  speed  Vj+i  is  set  to  zero  and  the  greatest 
speed  between  the  waypoints  vpk  is  calculated.  If  vpk  is  greater  than  the  specification  vcra,  then  vpk 
is  set  to  the  value  of  vcm;  ti,  t2,  and  taj+i  are  calculated;  and  the  interpolated  points  are  calculated.  If 
vpk  is  less  than  the  specification  vcra,  then  vpk  remains  unchanged,  and  t,  (again,  the  i  subscript 
represents  1  or  2  since  in  this  case  fi  =  t2)  is  calculated  with  taj+i,  and  the  interpolated  points  are 
also  calculated. 

If  the  next  waypoint  isn’t  the  final  point,  the  user  is  asked  to  input  the  latitude,  longitude,  and  (in 
the  event  the  platform  is  not  a  ground  platform)  the  elevation.  The  distance  between  the 
previous  waypoint  and  the  current  waypoint  Dj  is  calculated.  If  there  is  not  enough  distance  for 
the  platfonn  to  reach  the  cruising  velocity  vcru,  then  the  maximum  velocity  vmax  is  calculated; 
otherwise,  vmax  is  set  equal  to  vcru.  The  minimum  corresponding  arrival  time  tmin  is  calculated.  If 
there  is  not  enough  distance  for  the  platfonn  to  reach  zero  velocity,  the  minimum  velocity  vmm  is 
calculated,  and  the  loiter  time  is  set  to  zero  and  the  corresponding  value  of  tmax  calculated.  Then, 
Vmin,  tmax,  vmax,  and  tmin  are  displayed  as  guidance.  Otherwise,  vmin  is  set  equal  to  zero,  with  the 
values  tmin,  tmax,  vmin,  and  Vmax  displayed  as  guidance.  If  vmin  is  zero,  the  user  is  prompted  to 
input  the  loiter  time. 

If  the  loiter  time  is  zero,  then  final  data  for  the  platform  motion  at  the  waypoint  is  inputted,  either 
as  an  anival  time  taj+i  or  waypoint  speed  Vj+i.  The  deployment  module  calculates  the  maximum 
possible  speed  between  the  cunent  and  previous  waypoint  vtop  and  compares  it  to  the  cruise 
speed  vcru.  If  vtop  is  greater  than  vcru,  then  vtop  is  set  equal  to  vcru.  When  vtop  is  used,  the 
inflection  points  fi  and  t2  are  calculated,  and  the  departure  time  taj+i,  after  which  the  index  j  is 
increased  by  1,  and  the  index  are  tested  to  see  if  the  next  waypoint  is  the  final  waypoint. 
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4.2  Modifying  Input 

Next,  we  assume  that  a  string  of  waypoints  has  been  inputted  and  recorded  by  the  deployment 
module,  but  that  the  user  wishes  to  modify  them.  Four  processes  are  possible:  moving  a 
waypoint,  adding  a  waypoint,  removing  a  waypoint,  and  changing  either  the  arrival  time  taj+i, 
waypoint  velocity  Vj,  or  if  the  waypoint  velocity  Vj  is  zero,  the  loiter  time. 

4.2.1  Moving  a  Waypoint 

Figure  34  shows  the  sequence  for  moving  a  waypoint.  It  is  similar  to,  but  somewhat  more 
abbreviated  than,  the  process  for  inputting  a  waypoint.  After  starting  the  process,  the  user 
identifies  the  platfonn  and  waypoint  to  be  moved.  If  this  is  a  ground  platform,  the  user  inputs 
the  new  latitude  and  longitude  while  the  software  provides  the  elevation  data  from  the  DTED 
database  (22).  If  it  is  not  a  ground  platform,  the  user  provides  the  latitude,  longitude,  and 
elevation.  If  this  is  the  first  waypoint,  the  index  is  immediately  changed  because  the  waypoint 
speed  and  the  initial  time  are  zero  and  must  not  be  changed.  The  loiter  time  can  be  changed  by 
another  process.  Implicit  in  all  these  processes  is  the  assumption  that  every  platform  has  at  least 
two  distinct  waypoints.  The  case  where  a  platfonn  has  only  one  waypoint  will  be  treated  later. 

If  this  is  the  final  waypoint  for  the  platfonn,  the  distance  Dj  is  calculated,  and  if  there  is  enough 
distance  for  the  platform  to  stop  (2  a0  Dj  >  Vj  ),  the  final  speed  Vj+i  is  assigned  to  zero  and  the 
peak  speed  vpk  is  calculated.  If  the  platform  peak  speed  vpk  exceeds  the  cruising  speed  vcni,  then 
vpk  is  set  to  vcru,  and  the  anival  times  taj+i,  ti,  and  t2  (see  figure  29)  are  calculated  and  the 
interpolation  redone.  If  the  platform  is  not  greater  than  vcm,  then  the  single  inflection  time  t,  is 
calculated  (in  this  case,  fi  =  t2,  so  the  i  represents  the  subscripts  1  and  2,  since  there  is  only  one 
inflection  point)  along  with  the  arrival  time  taj+i  (see  figure  30),  after  which  the  interpolation  is 
redone.  Should  the  distance  Dj  not  be  sufficient  to  allow  the  platform  to  stop,  vmin  is  calculated 
and  the  speed  Vj+i  is  assigned  the  value  vmin;  the  arrival  time  taj+i  is  calculated,  along  with  the 
inflection  time  f  (in  this  case,  fi  =  t2,  so  the  i  represents  the  subscripts  1  and  2);  ti0iter  is  assigned 
to  0;  and  the  index  j  is  increased  by  one,  as  is  the  recording  index  (which  is  the  number  of 
waypoints  for  this  platform).  Then  data  is  input  for  the  next  waypoint  (latitude,  longitude,  and 
elevation  if  it  is  not  a  ground  platform,  and  latitude  and  longitude  with  elevation  read  from  a 
database  if  it  is  a  ground  platform),  and  the  distance  Dj  calculated.  The  distance  is  then  tested  to 
see  if  the  platfonn  has  enough  room  to  stop.  If  it  is  not,  the  user  is  prompted  for  another  point 
until  the  platform  has  distance  enough  to  stop. 

If  this  is  not  the  final  waypoint,  the  distance  Dj  is  calculated.  If  there  is  enough  distance  Dj  for 
the  speed  to  reach  vcru,  then  the  maximum  speed  of  the  platform  between  waypoints  vmax  is 
assigned  the  value  vcru.  If  there  is  not  enough  distance,  then  the  maximum  speed  the  distance 
allows  (equation  104)  is  assigned  to  vmax.  Next,  the  speed  for  this  waypoint  assigned  previously 
(Vj+i)  is  tested  to  see  if  it  is  greater  than  the  vmax.  If  it  is,  then  Vj+i  is  assigned  the  value  vmax  and 
if  it  is  not,  Vj+i  remains  unchanged.  Then,  the  distance  is  tested  to  see  if  the  platfonn  has  enough 
space  to  reach  zero  velocity.  If  it  does,  then  vm;n  is  assigned  the  value  of  0.  If  it  does  not,  then 
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Figure  34.  Overview  of  the  process  of  moving  a  waypoint. 


50 


2 

the  minimum  possible  velocity  vmin  is  assigned  the  lowest  speed  the  platfonn  can  reach  (vf  -  2  a0 
Dj)  “  by  equation  1 16.  Vj+i  is  next  tested  to  see  if  it  is  less  than  vmin.  If  it  is,  it  is  assigned  the 
value  of  vmin.  If  it  is  not,  it  remains  unchanged.  Next,  vtop  is  calculated  by  equation  122.  If  vtop 
is  greater  than  the  cruise  speed  vcra,  vtop  is  assigned  the  value  of  vcra.  If  it  is  not,  it  remains 
unchanged.  Next,  if  the  waypoint  velocity  Vj+i  is  nonzero,  the  loiter  time  is  set  to  zero.  If  Vj+i  is 
zero,  then  the  loiter  time  remains  unchanged.  Finally,  taj+i,  ti,  t2,  and  tdj+i  are  calculated,  and  the 
index  j  is  increased  to  check  the  next  waypoint. 

4.2.2  Adding  a  Waypoint 

The  flowchart  for  adding  a  waypoint  is  shown  in  figure  35.  Portions  of  the  figure  resemble  the 
process  of  moving  a  waypoint,  shown  in  figure  34.  After  identifying  the  platfonn  type  and  new 
index  number,  the  user  inputs  the  longitude,  latitude,  and  elevation  for  an  aerial  platform,  or  just 
the  latitude  and  longitude  for  a  ground  platform.  The  indices  for  the  waypoints  that  follow  the 
new  waypoint  in  time  are  incremented  by  one,  as  is  the  recording  index  that  records  the  number 
of  waypoints.  If  the  added  waypoint  is  intended  to  be  the  first  waypoint,  the  software  assigns  the 
waypoint  speed  and  arrival  time  as  zero.  The  user  is  prompted  for  the  loiter  time,  the  departure 
time  td,i  is  set  to  the  loiter  time,  and  the  index  is  increased  by  one.  The  next  waypoint  is  tested  to 
see  if  it  is  the  final  waypoint. 

If  the  added  waypoint  is  not  intended  to  be  the  first  waypoint,  it  is  tested  to  see  if  it  is  the  final 
waypoint.  Should  this  be  the  final  waypoint,  the  distance  Dj  is  calculated.  If  there  is  enough 
distance  for  the  platform  to  stop  (2  a0  Dj  >  Vj  ),  the  final  waypoint  velocity  Vj+i  is  assigned  to  zero 
and  the  peak  speed  vpk  is  calculated.  If  the  peak  speed  vpk  has  time  to  get  to  vcra,  then  vpk  is 
assigned  the  value  of  vcru,  and  the  arrival  times  taj+i,  ti,  and  t2  are  calculated  (see  figure  29)  and 
the  interpolation  redone.  Otherwise  (see  figure  30),  vpk  remains  unchanged  and  t,  is  calculated 
(the  i  subscript  represents  1  or  2,  since  in  this  case  ti  =  t2)  with  taj+i  and  the  interpolated  points 
are  recalculated.  If  there  is  not  enough  distance  for  the  platform  to  stop,  vmm  is  calculated  and 
the  speed  Vj+i  is  assigned  the  minimum  possible  value  for  speed  vmi„;  tmax  is  calculated  and 
assigned  to  the  arrival  time  taj+i;  hoiter  is  assigned  to  0;  and  the  index  j  is  increased  by  one,  as  is 
the  recording  index  (which  is  the  number  of  waypoints).  The  user  is  prompted  to  input  the 
latitude,  longitude,  and  elevation  for  an  aerial  platform,  or  only  latitude  and  longitude  for  a 
ground  platfonn.  The  distance  Dj  to  the  new  platform  is  calculated,  and  the  cycle  of  new 
platforms  is  inputted  until  it  is  able  to  stop. 

If  this  is  not  the  final  waypoint,  then  guidance  must  be  calculated  to  let  the  user  select  the  new 
waypoint  arrival  time  taj+i  or  waypoint  velocity  Vj+i.  This  is  first  accomplished  by  calculating 
the  distance  to  the  new  waypoint  Dj.  If  there  is  enough  distance  for  the  platform  to  reach  vcru, 

(2  a0  Dj  +  Vj  )  ,  then  vmax  is  assigned  the  value  of  vcru,  and  the  corresponding  tmjn  is  calculated. 

If  there  isn’t  enough  distance,  then  maximum  speed  vmax  is  calculated,  with  the  corresponding 

tmin- 
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Figure  35.  Overview  of  the  process  of  adding  a  waypoint. 
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Next,  the  guidance  values  for  vmin  and  tmax  are  calculated.  If  there  is  enough  distance  Dj  for  the 

2 

platform  to  stop  (2  a0  D,  -Vj  >  0),  then  vmin  is  assigned  a  value  of  zero  and  the  corresponding 
value  of  tmax  is  calculated.  The  formula  to  use  for  finding  tmax  depends  on  whether  there  is 
enough  distance  for  the  speed  to  reach  vcru.  If  there  is  [  vcru  <  (2  a0  Dj  +  vf  )  ],  then  tmax  is 

calculated  using  equation  106.  If  there  is  not,  tmax  is  calculated  with  equation  115,  followed  by 
the  display  of  vmin,  tmax,  vmax,  and  tmin  to  the  user  as  guidance,  who  inputs  the  loiter  time.  If  there 
is  not  enough  distance  Dj  for  the  platform  to  stop  (2  a0  Dj  -  Vj  <  0),  then  vmin  is  calculated  with 
equation  1 16,  and  the  corresponding  value  of  tmax  is  calculated  with  equation  1 17  and  displayed 
to  the  user  as  guidance  while  the  loiter  time  is  set  to  zero. 

If  the  loiter  time  is  not  zero,  the  waypoint  speed  Vj+i  is  set  to  zero.  The  highest  speed  possible  for 
the  platform  between  waypoints  vtop  is  calculated  with  equation  122,  and  if  it  exceeds  vcru,  vtop  is 
set  to  vcru,  and  taj+i  is  computed  with  equation  128.  If  the  loiter  time  is  zero,  the  user  is  given  the 
option  of  entering  the  arrival  time  ta,j+i  or  the  waypoint  velocity  Vj+i.  If  the  user  chooses  to  enter 
taj+i,  then  vtop  is  calculated  using  equation  125,  and  if  it  exceeds  vcra,  vtop  is  set  to  vcm,  and  Vj+i  is 
computed  using  equation  129.  If  the  user  chooses  to  enter  Vj+i,  then  vtop  is  calculated  with 
equation  122,  and  if  it  exceeds  vcru,  it  is  set  to  vcra,  and  taj+i  is  computed  with  equation  128.  Then 
ti  and  t2  are  calculated  with  equations  126  and  127,  taj+i  is  calculated,  and  the  waypoint  index  is 
increased  by  1 . 

If  the  index  is  for  the  final  waypoint,  then  it  is  treated  as  the  final  waypoint  as  outlined 
previously.  If  it  is  not,  it  and  all  the  remaining  waypoints  will  have  to  be  checked  for  self- 
consistency.  To  do  this,  the  distance  Dj  is  calculated,  and  from  this,  vmax  is  calculated  with 
equation  104.  If  vmax  exceeds  vcra,  vmax  is  set  to  vcru,  otherwise  it  retains  the  same  value.  The 
waypoint  speed,  Vj+i,  is  compared  to  vmax;  if  it  is  greater  than  vmax,  it  is  set  to  vmax.  If  it  is  not,  it 
retains  its  value.  Similarly,  if  2  a0  Dj  -  Vj  is  positive,  vmin  is  set  equal  to  zero,  or  the  square  root 
of  its  negative  if  2  a0  Dj  -  Vj  is  negative.  Then  Vj+i  is  compared  to  vmin  and  set  equal  to  vmin  if 
Vj+i  is  less  than  vmin-  A  new  variable,  vtop  is  computed  with  equation  122.  If  vtop  is  greater  than 
the  cruising  speed,  vcru,  vtop  is  set  equal  to  vcru.  If  Vj+i  is  greater  than  zero,  the  loiter  time  is  set  to 
zero.  Then  taj+i,  ti,  t2,  and  tdj+i  are  calculated,  and  the  waypoint  index  is  increased  by  1. 

4.2.3  Removing  a  Waypoint 

Removing  a  waypoint  is  a  shorter  procedure  and  is  illustrated  in  figure  36.  The  user  first 
identifies  the  index  and  type  of  platfonn  for  which  the  waypoint  must  be  removed.  After  reading 
the  acceleration  factor  a0  and  the  top  cruising  speed  vcra  from  a  database,  the  remaining  indices 
for  the  waypoints  for  this  platform  have  1  subtracted  from  them,  and  the  recording  index  has  1 
subtracted  from  it. 

If  the  waypoint  removed  was  the  first,  the  former  second  waypoint  becomes  the  new  first 
waypoint.  Then,  the  waypoint  velocity  and  arrival  time  are  set  to  zero.  The  user  is  prompted  for 
the  loiter  time,  and  the  index  is  increased  by  one.  The  next  waypoint  is  tested  to  see  if  it  is  the 
last  waypoint. 
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Figure  36.  Overview  of  the  process  of  removing  a  waypoint. 
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If  this  is  the  final  waypoint,  the  distance  Dj  is  calculated,  and  if  there  is  enough  distance  for  the 
platform  to  stop  (2  a0  Dj  >  Vj  ),  the  final  waypoint  velocity  Vj+i  is  assigned  to  zero  and  the  peak 
speed  vpk  is  calculated  with  equation  113.  If  vpk  is  greater  than  vcm,  then  vpk  is  assigned  the  value 
of  vcru,  the  arrival  times  taj+i,  ti,  and  to  are  calculated  (see  figure  29),  and  the  interpolation  is 
redone.  Otherwise  (see  figure  30),  vpk  remains  unchanged  and  t,  (the  i  subscript  represents  1  or 
2,  since  in  this  case  fi  =  t2)  and  taj+i  are  calculated  and  the  interpolated  points  recalculated.  If 
there  is  not  enough  distance  for  the  platfonn  to  stop,  the  speed  Vj+i  is  assigned  the  minimum 
possible  value  for  speed  vmin  after  calculating  vmin  and  the  arrival  time  taj+i,  while  ti0iter  is 
assigned  to  0  and  the  index  j  is  increased  by  one,  as  is  the  recording  index  (which  is  the  number 
of  waypoints).  The  user  is  prompted  to  input  the  latitude,  longitude,  and  elevation  for  an  aerial 
platform,  or  only  latitude  and  longitude  for  a  ground  platform.  The  distance  Dj  to  the  new 
platform  is  calculated  and  the  cycle  of  new  platforms  inputted  until  it  is  able  to  come  to  a  stop. 

If  it  is  not  the  final  waypoint,  it  and  all  the  remaining  waypoints  will  have  to  be  checked  for  self- 
consistency.  To  do  this,  the  distance  Dj  is  calculated,  and  from  this,  a  test  is  done  to  see  if  the 
maximum  speed  (2  a0  Dj  +  Vj  )  exceeds  vcra.  If  it  does,  vmax  is  set  to  vcra-  The  waypoint  speed 
Vj+i  is  compared  to  vmax;  if  it  is  greater  than  vmax,  it  is  set  to  vmax.  Similarly,  if  2  a0  Dj  -  Vj  is 
positive,  Vmin  is  set  equal  to  0,  or  if  2  a0  Dj  -  Vj  is  negative,  vmin  is  set  equal  to  the  square  root  of 
Vj'  -  2  a0  Dj.  Then  Vj+i  is  compared  to  vmin  and  set  equal  to  vmin  if  Vj+i  is  less  than  vmm.  A  new 
quantity,  vtop,  is  calculated,  and  if  it  exceeds  the  cruising  speed  vcru,  it  is  set  to  vcra-  If  vj+i  is  non¬ 
zero,  the  loiter  time  is  set  to  zero.  Then  taj+i,  ti,  t2  and  tdj+i  are  calculated,  and  the  waypoint 
index  is  increased  by  1 . 

4.2.4  Changing  Waypoint  Data 

The  final  process  involves  changing  the  waypoint  speed,  the  arrival  time  taj+i,  or  the  loiter  time 
of  a  single  waypoint.  Figure  37  summarizes  this.  The  user  first  identifies  the  index  and  type  of 
platform  for  which  the  waypoint  data  is  to  be  changed.  After  reading  the  acceleration  factor  a0 
and  the  top  cruising  speed  vcru  from  a  database,  the  waypoint  index  is  tested  to  see  if  it  is  the 
initial  point.  If  it  is,  the  process  prompts  the  user  to  input  a  new  loiter  time.  The  arrival  and 
departure  times  are  then  changed  for  the  remaining  waypoints  by  the  difference  of  the  old  and 
new  loiter  times  in  the  first  point.  It  is  unnecessary  to  change  the  other  velocities  because  the 
first  point  velocity  must  remain  zero.  The  process  then  stops.  Should  this  not  be  the  initial 
waypoint,  the  index  is  tested  to  see  if  it  is  the  final  waypoint.  If  it  is,  the  process  stops  since  the 
last  waypoint  must  have  a  waypoint  velocity  Vj+i  of  zero,  and  it  is  unnecessary  to  change  the 
arrival  time. 

If  the  index  is  not  for  the  final  waypoint,  the  distance  from  the  last  waypoint  Dj  is  recalculated. 

If  there  is  enough  distance  for  the  platform  to  reach  vcru,  (2  a0  Dj  +  Vj  )  ',  then  vmax  is  assigned 
the  value  of  vcni,  and  the  corresponding  tm;n  is  calculated  with  equation  101.  If  there  isn’t  enough 
distance,  then  maximum  speed  vmax  is  calculated  with  equation  104,  with  the  corresponding  tmm 
calculated  with  equation  105.  Next,  if  there  is  not  enough  distance  for  the  platform  to  reach  zero 
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Figure  37.  Overview  of  the  process  of  changing  the  waypoint  data. 
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velocity,  the  minimum  velocity  vmin  is  calculated  with  equation  1 16,  with  the  corresponding  tmax 
calculated  with  equation  1 17,  and  the  loiter  time  is  set  to  zero,  and  tmin,  tmax,  vmin,  and  vmax  are 
displayed  to  the  user  as  guidance.  Otherwise,  vmin  is  set  equal  to  zero,  and  the  maximum 
corresponding  arrival  time  tmax  is  calculated  with  either  equation  106  (if  the  maximum  speed  vmax 
is  vcru)  or  with  equation  115  (if  the  maximum  speed  vmax  is  vcra).  The  values  tmin.  tmax,  vmin,  and 
vmax  are  displayed  to  the  user  as  guidance.  If  vmin  is  zero,  the  user  is  prompted  to  input  the  loiter 
time. 

The  loiter  time  is  then  tested.  If  the  loiter  time  is  not  zero,  the  waypoint  speed  Vj+i  is  set  to  zero. 
The  highest  speed  possible  for  the  platform  between  waypoints  vtop  is  calculated  with  equation 
122,  and  if  it  exceeds  vcm,  it  is  set  to  vcra,  and  taj+i  is  computed  with  equation  128.  If  the  loiter 
time  is  zero,  the  user  is  given  the  option  of  entering  the  arrival  time  taj+i  or  the  waypoint  velocity, 
vj+i.  If  the  user  chooses  to  enter  taj+i,  then  vtop  is  calculated  with  equation  125,  and  if  it  exceeds 
Vcm,  it  is  set  to  vcm,  and  Vj+i  is  computed  with  equation  129.  If  the  user  chooses  to  enter  Vj+i,  then 
again  vtop  is  calculated  with  equation  122,  and  if  it  exceeds  vcm,  it  is  set  to  vcra,  and  taj+i  is 
computed  with  equation  128.  Then  ti,  h,  and  tdj+i  are  calculated  using  equations  126  and  127, 
and  the  waypoint  index  is  increased  by  1 . 

The  index  is  again  tested  to  see  if  it  corresponds  to  the  final  waypoint.  If  this  is  the  final 
waypoint,  the  distance  Dj  is  calculated,  and  if  there  is  enough  distance  for  the  platform  to  stop  (2 
a0  Dj  >  Vj  ),  the  final  waypoint  velocity  Vj+i  is  assigned  to  zero  and  the  peak  speed  vpk  is 
calculated  with  equation  113.  If  the  peak  speed  vpk  exceeds  vcru,  then  vpk  is  assigned  the  value  of 
Vcm,  and  the  arrival  times  taj+i,  ti,  and  t2  are  calculated  with  equations  75,  107,  and  108  (see 
figure  29)  and  the  interpolation  is  redone.  Otherwise  (see  figure  30),  vpk  remains  unchanged  and 
h  (the  i  subscript  represents  1  or  2,  since  in  this  case  fi  =  to)  and  taj+i  are  calculated  with 
equations  110  and  111,  and  the  interpolated  points  are  recalculated.  If  there  is  not  enough 
distance  for  the  platform  to  stop,  the  speed  vmin  is  calculated  with  equation  1 16  and  the  speed  Vj+i 
is  assigned  the  value  for  vmin,  the  arrival  time  taj+i  (which  equals  tmax)  is  calculated  with  equation 
1 17,  ti0iter  is  assigned  to  0,  and  the  index  j  is  increased  by  one,  as  is  the  recording  index  (which  is 
the  number  of  waypoints).  The  user  is  prompted  to  input  the  latitude,  longitude,  and  elevation 
for  an  aerial  platfonn,  or  only  latitude  and  longitude  for  a  ground  platform.  The  distance  Dj  to 
the  new  platform  is  calculated  and  the  cycle  of  new  platforms  inputted  until  the  platform  is  able 
to  come  to  a  stop. 

If  it  is  not  the  final  waypoint,  it  and  all  the  remaining  waypoints  will  have  to  be  checked  for  self- 
consistency.  To  do  this,  the  distance  Dj  is  calculated,  and  from  this,  a  test  is  done  to  see  if  the 
maximum  speed  (2  a0  Dj  +  Vj  )  can  exceed  vcra-  If  it  can,  vmax  is  set  to  vcra;  otherwise,  vmax  is 
calculated  with  equation  104.  The  waypoint  speed,  Vj+i,  is  compared  to  vmax — if  it  is  greater  than 
Vmax,  it  is  set  to  vmax.  Otherwise  Vj+i  remains  unchanged.  Similarly,  if  2  a0  Dj  -  vf  is  positive, 
vmin  is  set  equal  to  0,  or  (vj~  -  2  a0  Dj)  "  if  2  a0  Dj  -  vf  is  negative,  and  the  loiter  time  is  set  to 
zero.  Then  Vj+i  is  compared  to  vmin  and  set  equal  to  vmin  if  Vj+i  is  less  than  vmin-  Otherwise,  Vj+i 
remains  unchanged.  A  new  variable,  vtop  is  calculated  with  equation  122.  If  it  exceeds  vcru,  it  is 
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set  to  vcru.  Otherwise,  it  remains  unchanged.  If  Vj+i  is  nonzero,  the  loiter  time  is  set  to  zero.  Then 
taj+i,  ti,  t2,  and  tdj+i  are  calculated  with  equations  126-128,  the  waypoint  index  is  increased  by  1, 
and  the  departure  time  taj+i  is  calculated. 

4.2.5  If  a  Platform  Has  Only  One  Waypoint 

We  now  consider  the  input  process  and  the  four  modification  processes  (moving  a  waypoint, 
adding  a  waypoint,  removing  a  waypoint,  and  reassigning  values  to  the  waypoint)  when  a 
platform  has  only  one  waypoint  assigned  to  it.  When  a  platform  has  only  one  waypoint,  the 
waypoint  velocity  (vi)  is  assigned  to  zero.  The  arrival  time  ta,i  is  also  assigned  to  zero.  The 
departure  time  td,i  is  equal  to  the  loiter  time,  which  equals  the  mission  time.  If  a  platform  has 
only  one  waypoint,  we  assume  that  it  does  not  move  but  remains  stationary  for  the  entire  run.  As 
a  result,  the  interpolation  process  for  this  platfonn  need  only  copy  the  initial  waypoint  latitude, 
longitude,  and  elevation  for  every  snapshot  time  to  the  snapshot  database. 

Inputting  data  for  a  single  waypoint  platfonn  is  relatively  straightforward  and  is  set  out  in  figure 
38.  The  user  identifies  the  type  of  platform  that  requires  a  single  waypoint,  and  the  default  index 
(as  well  as  the  recording  index)  is  set  to  1 .  The  arrival  time  ta,i  and  the  waypoint  velocity  vi  are 
set  to  zero.  The  loiter  time  becomes  the  mission  time,  so  that  the  departure  time  td,i  is  set  equal 
to  the  loiter  time.  If  the  platform  is  a  ground  platform,  the  user  inputs  the  latitude  and  longitude, 
while  the  elevation  is  read  from  a  DTED  database.  If  it  is  not  a  ground  platform,  the  user  also 
inputs  the  elevation.  The  interpolation  is  done  by  copying  the  latitude,  longitude,  and  elevation 
database  to  the  snapshot  database  for  this  platform. 


Figure  38.  The  process  of  inputting  data  for  a  platform  with  one  waypoint. 


Moving  a  waypoint  (when  a  platfonn  has  a  single  point)  only  requires  a  reassignment  of  a  new 
latitude  and  longitude.  Once  the  user  identifies  the  platform,  the  index  is  assumed  to  be  1 .  If  it 
is  a  ground  platform,  the  user  inputs  the  new  latitude  and  longitude.  If  it  is  not,  the  user  inputs 
the  new  latitude,  longitude,  and  elevation.  Since  that  is  all  that  is  required,  the  process  ends. 
This  brief  process  is  outlined  in  figure  39. 
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Figure  39.  The  process  of  inputting  data  for  changing  the  location  of  a  single  waypoint  platform. 


Adding  a  point  is  a  bit  more  involved.  After  identifying  the  platform  on  which  to  add  a 
waypoint,  the  user  identifies  whether  the  added  point  comes  before  the  previously  defined  point 
in  time  (the  new  point  has  an  index  of  j  =  1)  or  after  it  (the  new  point  has  an  index  of  j  =  2).  The 
platform  acceleration  a0  and  cruising  speed  vcm  are  read  from  a  database,  and  the  recording  index 
is  set  to  2.  If  it  is  a  ground  platform,  the  user  inputs  the  new  latitude  and  longitude.  If  it  is  not, 
the  user  inputs  the  new  latitude,  longitude,  and  elevation.  The  two  points  then  (if  necessary) 
have  their  indices  modified  so  that  the  first  point  in  time  has  an  index  of  1  and  the  second 
waypoint  in  time  has  an  index  of  2.  Since  we  now  have  two  distinct  waypoints,  the  process 
continues  as  outlined  in  figure  35.  The  entire  process  is  laid  out  in  figure  40.  The  verbiage 
describing  the  remaining  part  of  the  process  is  identical  to  the  verbiage  used  to  describe  figure 
35.  For  the  sake  of  brevity,  it  will  not  be  repeated  here. 

Removing  a  data  point  for  a  platform  with  only  one  data  point  is  tantamount  to  deleting  the 
platform.  So  instead  of  outlining  a  separate  process  to  remove  the  only  waypoint  to  a  platform, 
the  process  will  be  to  simply  delete  the  platform.  Likewise,  reassigning  waypoint  data  for  a 
platform  with  only  one  waypoint  is  unnecessary  since  there  is  only  one  waypoint.  The  single 
waypoint  must  have  the  waypoint  velocity  Vi  and  arrival  time  ta,i  equal  to  zero,  so  the  data  can 
not  be  changed.  Consequently,  this  process  will  not  be  included  in  the  deployment  module. 

4.3  The  Interpolation  Process 

Once  the  latitudes,  longitudes,  and  elevations  have  been  established  for  the  set  of  waypoints  for  a 
platform  (with  their  corresponding  arrival  and  departure  times  either  assigned  directly  by  the 
user  or  calculated  from  assigned  waypoint  velocities  and  loiter  times),  then  it  is  possible  to 
calculate  the  platform’s  latitude,  longitude,  and  elevation  for  a  selected  snapshot  time  from  the 
interpolation  time  f.  The  process  is  to  first  calculate  the  interpolated  distance  D;  from  the  most 
recent  waypoint  the  platform  encountered  to  the  platfonn’s  position  at  f.  The  interpolated 
distance  D;  makes  it  possible  to  obtain  the  interpolated  elevation  Ei.  The  interpolated  angle  A9j 
that  the  interpolated  location  of  the  platform  at  time  t;  makes  with  the  location  of  the  last 
waypoint  at  the  Earth’s  center  then  makes  it  possible  to  calculate  the  interpolated  longitude  9i 
and  latitude  cj)j  after  making  a  coordinate  transformation  to  simplify  the  use  to  the  quantity  A9j. 
Each  of  these  processes  will  be  explained  in  detail. 
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Figure  40.  The  process  of  adding  a  waypoint  to  a  platform  that  has  one 
waypoint. 
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4.3.1  Calculate  the  Interpolated  Distance 

Assuming  that  we  know  the  interpolated  time  fi,  the  first  task  is  to  see  if  any  of  the  waypoint 
loiter  periods  (the  time  between  the  arrival  time  tj,a  and  departure  time  tdj)  bound  the  interpolated 
time  f.  Put  another  way,  is  taj  <  h  <  tdj?  If  it  is,  then  the  interpolated  latitude,  longitude,  and 
elevation  for  the  platform  are  set  equal  to  the  latitude,  longitude,  and  elevation  of  the  waypoint. 
This  ends  the  interpolation  process  for  the  interpolated  time  t,.  (Note  that  the  last  waypoint  loiter 
period  is  the  span  between  the  final  waypoint  arrival  time,  tN  a  and  the  ending  time  of  the  run.) 

If  the  interpolated  time  t;  does  not  occur  during  the  waypoint  loiter  time,  then  it  occurs  during  the 
transit  time  between  points.  The  interpolation  process  then  becomes  more  involved.  To 
calculate  the  interpolated  distance  D,,  it  is  first  necessary  to  compare  it  with  the  times  that  mark 
the  inflection  points  fi  and  t2  of  the  platform  profile  of  figure  24.  As  a  review,  we  recount  how 
the  inflection  points  are  calculated.  Step  1  is  to  calculate  vtop  from  equation  125,  restated  here  as 
equation  130. 

Vtop  =  Vj  +  a0  (taj+1  -  tdj  )  -  [  Vj  a0  (taj+1  -  tdj  )  -  a0  Dj  +  (  1  /  2  )  a02  ( taJ+i  -  tdj  )2  ] 1/2 .  (130) 

If  vtop  exceeds  vcra,  then  vtop  is  set  equal  to  vcru.  Otherwise,  equation  130  defines  vcru.  Equations 

126  and  127  (repeated  here  as  equations  131  and  132)  are  then  used  to  establish  the  inflection 
times  ti  and  t2. 

ti  =  tdj  +  (  vtop-  Vj)/a0.  (131) 

t2  =  taj+1  -  (  Vtop  -  Vj+1  )  /  a0 .  (1 32) 

To  calculate  the  distance  D;  from  the  most  recently  departed  waypoint,  it  is  necessary  to 
established  the  relationship  between  the  departure  time  tdj,  the  inflection  times  fi  and  t2,  the 
arrival  time  at  the  next  waypoint  taj+i,  and  the  interpolation  time  t,.  The  next  three  equations 
(133-135)  are  taken  from  equations  64,  67,  and  70. 

If  tdj  <  h  <  t1?  D,  =  a0  ( tf  -  tdj2 )/  2  -  (  a0  tdj  -  v, )  ( t,  -  tdj ) .  (133) 

If  ti  <  t;  <  t2,  Dj  =  a0  ( ti2  -  tdj2  )/  2  -  (  a0  tdj  -  vj )  (  fi  -  tdj  )  +  vcru  (t,  -  fi  )  .  (134) 

If  t2  <  h  <  taj+i,  D;  =  a0  ( ti2  -  tdj2 )/  2  -  (  a0  tdj  vj )  ( t,  tdj )  +  vcru  (t2  -  ti )  + 

Vcru  ( t;  -  t2  )  -  a0  ( t;2  - 122  )  /  2  +  a0  t2  ( t;  - 12)  .  (135) 

Now  that  we  have  the  interpolated  distance  Dj,  we  can  calculate  the  interpolated  elevation,  E,. 

4.3.2  Calculate  the  Interpolated  Elevation 

If  the  interpolated  time  t;  does  not  occur  during  the  waypoint  loiter  time,  it  is  possible  to 
calculate  the  interpolated  elevation  E,  as  a  result  of  one  of  two  circumstances.  First,  the  two 
waypoints  bounding  the  platform  in  time  have  the  same  elevation.  In  that  case,  Ej  =  Ej  =  Ej+i. 
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But  should  Ej  not  equal  Ej+i,  the  interpolated  elevation  must  be  found  using  a  two-step  process. 
First  reconsider  equation  97  by  eliminating  then  solving  for  the  elevation  constant  k: 

k  =  (Ej+i  -  Ej)  /  Dj  .  (136) 

Note  that  this  calculation  is  made  with  the  distance  between  the  waypoints  Dj,  not  the 
interpolated  distance  D;.  To  obtain  the  interpolated  elevation  Ei,  we  use  equation  97  again  but 
substitute  the  interpolated  distance  D,  (calculated  in  the  previous  section)  for  the  interwaypoint 
distance  Dj,  and  the  interpolated  elevation  E;  for  the  arrival  waypoint  elevation  Ej+i.  Rearranging 
the  terms  we  find 


k  D;  +  Ej  =  E; .  (137) 

With  the  interpolated  distance  D,  and  the  interpolated  elevation  Ej,  we  may  next  find  the 
interpolated  angle  between  the  line  connecting  the  most  recently  encountered  waypoint  and  the 
center  of  the  Earth,  and  the  line  connecting  the  platform’s  position  at  the  interpolated  time  t,  and 
the  center  of  the  Earth. 


4.3.3  Calculate  the  Interpolated  Angle  A0j 

Again,  we  assume  that  the  interpolated  time  f  does  not  occur  during  the  waypoint  loiter  time. 
And  again,  there  are  two  circumstances  we  need  to  consider  in  the  calculation  of  the  interpolated 
angle  A9;:  the  elevations  of  the  waypoints  bounding  the  interpolated  point  are  either  equal  or 
different. 


First,  if  the  elevations  of  the  waypoints  bounding  the  interpolated  point  in  time  are  equal,  we 
make  use  of  equation  99,  substituting  the  interpolated  angle  A9j  for  the  angle  the  waypoints  make 
at  the  Earth’s  center  A0,  and  the  distance  to  the  interpolated  point  Dj  for  the  distance  between  the 
waypoint  Dj.  When  we  solve  for  A0;, 

D;  /  (Rave  +  Ej)  =  A0; .  (138) 

Second,  if  the  elevations  of  the  waypoints  bounding  the  interpolate  point  in  time  are  not  equal, 
we  use  equation  94  instead.  Making  use  of  the  interpolated  elevation  Ei  and  the  constant  k,  and 
substituting  the  interpolated  elevation  Ej  for  the  waypoint  elevation  Ej+i  and  solving  for  A0i,  we 
find 


A0i  =  [(1  -  k2)1/2  /  k]  Ln  [(Rave  +  Ei)  /  (Rave  +  Ej  )]  .  (139) 

Ln  denotes  the  natural  log  function.  Note  that  if  the  platform  is  descending,  the  elevation 
constant  k  will  be  negative  (see  equation  136).  If  the  platform  is  descending,  then  Ej  will  be 
greater  than  Ei,  causing  the  natural  log  of  the  expression  to  be  negative.  We  assume  that  we  take 
the  positive  root  for  [(1  -  k  )  ,  resulting  in  a  positive  interpolation  angle  A0j.  Should  the 
platform  be  ascending,  k  and  the  natural  log  will  be  positive,  again  resulting  in  a  positive 
interpolation  angle  A9j.  Note  that  when  the  elevation  does  not  change  (the  situation  described  by 
equation  138),  the  interpolation  angle  A0,  is  positive  because  the  distance  interpolated  D,,  the 
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Earth’s  average  radius  Rave,  and  the  waypoint  elevation  Ej  are  always  positive.  The  fact  that  the 
interpolation  angle  AO,  is  always  positive  will  play  a  vital  role  in  determining  the  interpolated 
longitude  and  latitude. 

4.3.4  Calculate  the  Interpolated  Latitude  and  Longitude 

To  calculate  the  interpolated  longitude  and  latitude  involves  doing  a  coordination  transformation 
such  that  the  two  waypoints  bounding  the  interpolated  point  in  time  will  be  on  the  transformed 
equator.  Hence,  the  interpolated  point  will  also  be  on  the  transformed  equator,  which  will  be 
calculated  by  using  the  interpolated  angle  A0i.  By  reversing  the  process,  we  can  obtain  the 
interpolated  longitude  and  latitude. 

The  transformation  will  be  done  as  a  two-part  process.  The  first  part  involves  keeping  the 
latitude  constant  but  changing  the  longitude.  This  process  is  illustrated  in  figure  41 .  The  new  x 
axis  (denoted  x’  in  the  figure)  begins  at  the  Earth’s  center  and  intersects  the  point  on  the  equator 
where  the  line  of  constant  latitude  0o  crosses  the  equator.  The  y’  axis  begins  at  the  Earth’s 
center,  is  in  the  equatorial  plane,  and  is  90°  from  the  x’  axis.  The  axis  is  occulted  by  the  Earth  in 
the  figure,  although  the  original  y  axis  (also  in  the  equatorial  plane  and  beginning  at  the  Earth’s 
center)  is  not  occulted  and  is  plainly  shown  in  the  figure.  The  z  and  z’  axis  are  the  same  and  start 
at  the  Earth’s  center  and  progress  through  the  North  Pole. 

Recall  that  the  parameters  (t)max  and  0o  describe  the  equation  of  a  great  circle  and  are  derived  from 
the  latitudes  and  longitudes  of  the  two  bounding  waypoints  in  equations  10  and  11,  restated  here 
as  equations  140  and  141. 

0o  =  arc  tan  [(tan  (j)j  cos  0j+i  -  tan  4>j+i  cos  0j)  /  (tan  (j)j+i  sin  0j  -  tan  (j)j  sin  0j+i)]  .  (140) 

<j)max  =  arc  tan  [tan  4>j+i  /  cos  (0J+i  -  Oo)]  =  arc  tan  [tan  cj>j  /  cos  (0j  -  0o)]  .  (141) 

Recall,  too,  that  this  assumes  that  the  argument  for  arc  tan  is  finite  in  equation  140,  and  that  we 
take  the  positive  value  for  c[)max  for  the  value  of  0o  calculated  in  equation  140.  Additionally,  the 
point  (cj)max,Oo)  in  the  nonprimed  coordinate  system  is  the  point  of  highest  latitude  on  the  great 
circle.  Since  we  want  the  line  of  constant  latitude  that  goes  through  this  point  to  be  the  prime 
meridian  in  the  primed  system,  the  transformed  coordinate  will  be  (cjwkO0). 

So,  the  effect  of  the  first  part  of  the  transfonnation  will  be  simply  to  change  the  longitude  by  the 
quantity  0o  while  allowing  the  latitude  to  remain  unchanged.  A  way  to  visualize  this 
transformation  is  to  remove  the  Earth  and  look  along  the  z  axis  in  the  -z  direction.  Figure  42 
illustrates  the  transformation.  From  the  figure,  it  is  easy  to  obtain  the  transformation 
mathematically  for  not  only  the  interpolated  points,  but  for  any  point  that  includes  the  waypoints 
bounding  the  interpolated  point: 
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Figure  42.  The  first  step  in  transforming  the 
interpolated  point  to  a  new 
coordinate  system  with  the  Earth 
removed  viewing  in  the  -z 
direction. 
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0i’  =  0i-0o. 

(147) 

We  must  make  sure  that  the  three  transformed  longitudes  0j  ,  0j+i  ,  and  0i’  are  in  the  range  -180° 
to  +180°.  If  any  of  them  are  less  than  -180°,  we  add  360°.  If  any  of  them  exceed  +180°,  we 
subtract  360°. 

The  second  part  of  the  transformation  is  more  complicated.  It  involves  defining  a  set  of  double- 
primed  coordinates  such  that  the  primed  and  double -primed  y  axis  are  the  same.  Furthermore, 
the  x”  axis  begins  at  the  center  of  the  sphere  and  goes  through  the  point  ((|>max,0o)  in  the  original 
system  and  (())max,0)  in  the  single-primed  system.  This  line  intersects  the  double-primed  equator 
and  the  prime  meridian.  The  transformation  is  illustrated  in  figure  43. 


Figure  43.  The  second  step  in  transforming  the  interpolated  point  to  a  new  coordinate  system:  the  y’ 
axis  acts  as  the  axel  so  that  y’  =  y”. 
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To  obtain  the  quantitative  relationship  between  the  primed  and  double-primed  coordinate 
systems,  we  view  along  the  y’/y”  axis  in  the  positive  y’/y”  direction  with  the  Earth  removed  for 
clarity.  Figure  44  shows  this  situation.  The  fact  that  the  waypoints  and  the  interpolated  points 
he  on  the  equator  in  the  double-primed  coordinate  system  simplifies  the  transformation.  The 
next  few  equations  will  pertain  to  any  point  on  the  double-primed  equator.  As  the  figure  shows, 
any  point  on  the  equator  has  an  x”  coordinate  and  a  y”  coordinate.  Because  the  point  is  not  the 
transformed  double-primed  equator,  it  has  no  z”  coordinate.  Note  that  y’  =  y”,  so  that 
relationship  is  trivial.  Concentrating  on  the  relationship  between  the  x”  coordinate  and  the  z’ 
and  x’  coordinate,  we  see  from  the  figure  that 


z’  =  x”  sin  cj)max 

(148) 

x’  =  x”  cos  (jw . 

(149) 

Figure  44.  The  second  step  in  transforming  the  interpolated 
point  to  a  new  coordinate  system  with  the  Earth 
removed  viewing  in  the  +y’  direction. 

Making  use  of  equations  1-3  for  the  relationship  between  spherical  and  rectangular  coordinates, 
and  using  the  Ravc,  the  Earth’s  average  radius,  and  the  elevation  E,  we  may  rewrite  equations  148 
and  149  for  any  point  on  the  great  circle  as 

(  Rave  +  E  )  sin  4>’  =  (  Rave  +  E  )  cos  cj)”  cos  0”  sin  (j)max  (150) 

and 

(  Rave  +  E  )  COS  cj)’  COS  0’  =  (  Rave  +  E  )  COS  (j)”  COS  0”  COS  (j)max  ,  (151) 
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and  the  fact  that  the  points  on  the  double-primed  equator  have  (j>”  =  0  and  we  may  divide  by  the 
common  factor  (  Rave  +  E  )  renders  the  relationships  between  the  primed  and  double-primed 
latitudes  and  longitudes  as 


and 


sin  (j)’  =  cos  0”  sin  (j)max 

(152) 

COS  (j)’  cos  0’  =  cos  0”  cos  c()max . 

(153) 

Substituting  the  relationships  between  the  original  and  the  double-primed  coordinates  as  outlined 
in  equations  142  to  147,  and  solving  for  the  double-primed  longitude  0”,  we  find  for  the  jth 
waypoint,  the  j+ 1th  waypoint,  and  the  interpolated  point  such  that 


0j”  =  arc  cos  (  sin  4>j  /  sin  (j)max  )  .  (154) 

0j”  =  arc  cos  [  cos  cj)j  cos  (  0j  -  0o)  /  cos  (t)max  ]  ■  (155) 

0j+i”  =  arc  cos  (  sin  cj)j+i  /  sin  <|)max  )  .  (156) 

0j+i”  =  arc  cos  [  cos  cj)j+i  cos  (  0j+i  -  0o)  /  cos  cj)max  ]  .  (157) 

0i”  =  arc  cos  (  sin  (j);  /  sin  cj)max  )  .  (158) 

0”  =  arc  cos  [  cos  cj)i  cos  (  0;  -  0o)  /  cos  c))max  ]  .  (159) 


The  next  relationship  needed  for  the  transfonnation  takes  advantage  of  the  fact  that  the  primed 
and  double-primed  coordinate  systems  share  the  same  y  axis.  So  in  general,  y’  =  y”.  Using  the 
definition  of  equation  2,  we  may  establish  that 


y’  =  (  Rave  +  E  )  cos  c()'  sin  0’  =  y”  =  (  Rave  +  E  )  cos  (j)”  sin  0”.  (160) 


For  the  case  where  the  point  is  on  the  double-primed  equator,  we  may  take  (j)”  as  zero.  After 
dividing  by  (  Rave  +  E  ),  we  see  for  a  point  on  the  double-primed  equator 


cos  (j)' sin  0' =  sin  0”.  (161) 

Again,  substituting  the  values  for  the  jth  and  j+lth  waypoints  and  the  interpolated  point  from 
equations  142  to  147,  and  solving  for  the  longitude  in  the  double-primed  coordinate  system,  we 
find  that 


0”j  =  arc  sin  [cos  <j)j  sin  (  0,  -  0O  )  ]  .  (162) 

0”j+i  =  arc  sin  [cos  (j)J+i  sin  (  0  j+i  —  0O  )  ]  •  (163) 

0”i  =  arc  sin  [cos  (j);  sin  (  0;  -  0O  )  ]  .  (164) 

Note  that  equations  162-164  use  the  inverse  sine,  or  arc  sin,  function  to  obtain  the  longitude  in 
the  double-primed  coordinate  system  0”,  while  equations  154-159  use  the  inverse  cosine,  or  arc 
cos,  function  to  get  0”.  At  first  glance,  these  two  functions  appear  redundant.  But  the  inverse 
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sine  and  cosine  functions  deliver  not  one  but  two  answers  in  the  range  -180°  to  +180°.  This  is 
illustrated  in  figure  45  (28).  As  an  example,  dots  are  placed  on  the  sine  and  cosine  functions 
equal  to  0.5,  which  correspond  to  an  angles  of  30°  and  150°  for  sine,  and  60°  and  -60°  for  cosine. 


X  (Degrees) 


Figure  45.  The  dots  on  the  sine  and  cosine  graph  at  0.5  show  that  there  are  two  distinct  angles  for  any  one 
value  of  sine  and  cosine  in  the  interval  -180°  to  +180°. 
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To  sort  out  which  of  the  two  angles  that  each  of  the  inverse  functions  produce  is  the  correct 
angle,  we  must  first  examine  the  sign  of  the  functions  in  each  of  the  four  quadrants.  This  is 
shown  in  figure  46.  In  the  case  where  the  angle  9  is  between  0°  and  90°  (quadrant  I),  then  cos  9 
and  sin  9  are  both  positive.  When  9  is  between  99°  and  189°  (quadrant  II),  then  sin  9  is  positive 
again  but  cos  9  is  negative.  If  9  is  between  9°  and  -99°  (quadrant  IV),  then  cos  9  is  positive  but 
sin  9  is  negative.  Finally,  if  9  is  between  -99°  and  -189°  (quadrant  III),  then  both  cos  9  and  sin 
9  are  negative. 


Quadrant  II 

90°  <  0E  <  1 80° 

Cos  0E  <  0 

Sin  0E  >  0 

NN^n" 

Quadrant  I 

0°  <  6l  <  90° 

Cos  >  0 

Sin  6l  >  0 

n'x^iw/ 

Quadrant  III 

Quadrant  IV 

-90°  >  0m>- 180° 

o 

o 

V 

< 

V 

CD 

o 

o 

Cos  0E1  <  0 

Cos  tfIV  >  0 

Sin  0El<  0 

Sin  01V<  0 

Figure  46.  For  an  angle  9  ranging  from  -180°  to  +180°,  the  sine  and  cosine  functions  take 
on  either  a  positive  or  negative  sign  in  each  90°  quadrant. 


Next,  we  must  consider  which  values  of  9  are  returned  by  programming  languages  for  the 
functions  arc  cos  and  arc  sin.  As  is  the  case  for  Java  (29),  C++  (30,  31),  and  FORTRAN  (the 
Formula  Translator)  (32),  9  =  arc  cos  x  returns  a  value  for  9  that  is  in  the  first  and  second 
quadrants,  that  is  9  >  9  >  189°,  while  9  =  arc  sin  x  returns  a  value  for  9  that  is  in  the  first  and 
fourth  quadrants,  that  is  -99°  <  9  <  +99°.  Yet  9  could  be  in  any  of  the  four  quadrants  of  figure 
46.  We  see  that  for  all  angles  of  9  for  which  sin  9  is  positive  (those  in  quadrants  I  and  II),  the  arc 
cos  function  returns  the  value  of  9  correctly.  For  the  other  two  quadrants  (III  and  IV),  we  see 
that  the  sine  function  is  negative,  but  the  arc  cos  function  returns  a  value  of  9  that  is  opposite  in 
sign  to  the  values  nominally  returned  using  programming  languages.  Therefore,  to  make  sure  the 
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angle  9  returned  is  the  correct  case,  we  leave  the  sign  for  9  as  is  for  the  function  arc  cos  x  if  the 
sin  function  is  positive  and  change  the  sign  for  9  =  arc  cos  X  if  the  sin  function  is  negative. 
Applying  equations  154-159  and  162-164,  we  find  that 

if  cos  4>j  sin  (  9j  -  9o  )  >  9  , 

9j”  =  arc  cos  (  sin  4>j  /  sin  4>max  )  , 

and 

if  cos  (j)j  sin  (  9j  -  90  )  <  9  , 

9j”  =  -arc  cos  (  sin  (j)j  /  sin  cj)max  )  .  (165) 

Likewise, 

if  cos  (j)j+i  sin  (  9J+i  -  9o  )  >  9 
9j+i”  =  arc  cos  (  sin  cj)j+i  /  sin  ())max  )  , 
and 

if  cos  cj)j+i  sin  (  9j+i  -  90  )  <  9 

9j+i”  =  -arc  cos  (  sin  (j)j+i  /  sin  (j)max  )  .  (166) 

And  finally, 

if  cos  (|>i  sin  (  9;  -  90  )  >  9 
9;”  =  arc  cos  (  sin  (j);  /  sin  4>max  )  , 
and 


if  cos  4>i  sin  (  9i  -  9o  )  <  9 

9”  =  -arc  cos  (  sin  (j)i /  sin  c()max  )  .  (167) 

We  now  apply  equations  138  and  139  to  find  the  interpolated  angle  9”.  They  are  repeated  here 
as  equations  168  and  169. 

A9i  =  Di  /  (Rave  +  Ej)  .  (168) 

A9i  =  [(1  -  k2)1/2  /  k]  Ln  [(Rave  +  EO  /  (Rave  +  E,  )]  .  (169) 

Equation  168  is  used  if  the  elevations  at  the  waypoints  Ej  and  Ej+i  are  equal.  If  they  are  not,  then 
equation  169  is  used,  where  the  value  for  k  is  defined  by  equation  136,  repeated  here  as  equation 
179: 


k  =  (Ej+i  -  Ej)  /  Dj . 


(179) 


79 


Notice  that  it  is  Dj  in  equation  170,  the  distance  between  the  waypoints,  not  to  be  confused  with 
the  interpolated  distance  Di,  the  distance  from  the  jth  waypoint  to  the  platform  at  the  interpolated 
time  h.  One  more  item  needed  is  the  angular  distance  between  the  waypoints.  That  is,  the  angle 
between  two  lines:  one  that  goes  from  waypoint  j  to  the  Earth’s  center,  and  the  other  that  goes 
between  the  j+ 1th  waypoint  to  the  Earth’s  center.  This  angle,  A9,  was  calculated  in  equation  80 
and  is  repeated  here  as  equation  171. 

2 

A0  =  2  arc  sin  {  [(  1  /  2  )  [(cos  4>j+i  cos  9j+i  -  cos4>j  cos  0j)  + 

2  2  1/2 

(cos  4>j+i  sin  0j+i  -  cos  4>j  sin  0j)“  +  (sin  c[)j+i  -  sin  (j)j)“]  }  .  (171) 

(c()j,0j)  and  (4>j+i,0j+1)  are  the  longitude  and  latitude  of  the  jth  and  j+lth  waypoints  in  the  original 
coordinate  system.  The  interpolated  angle  in  the  double-primed  coordinate  system  9”  is  the 
fraction  of  the  angular  distance  A9  between  the  waypoints  that  the  platform  reached  at  the 
interpolated  time  t,.  Expressed  quantitatively  using  the  terms  derived  so  far,  it  is  the  fraction 
A9;  /  A9.  The  radial  distance  covered  between  the  waypoints  by  the  platfonn  is  the  fraction  of 
the  distance  covered  times  the  difference  between  the  coordinates  in  the  double-primed 
coordinate  system,  or  (  Qj+i”-  6j”  )  (  A0;  /  A0  ).  To  obtain  the  double-primed  interpolated 
longitude,  we  need  only  add  this  expression  to  the  longitude  of  the  jth  waypoint. 

0i”  =  0j”  +  (  0j+r-  0j”  )  (  A0i  /  A0  )  .  (172) 

Note  that  the  quantities  A0,  and  A9  are  always  positive.  But  depending  on  the  direction  of  travel, 
0j”  may  be  greater  than  or  less  than  0j+i”.  To  make  sure  that  0]  is  between  0(  and  9j+i  ,  it  is 
necessary  to  use  the  full  expression  in  equation  172  rather  than  simply  add  A0,  to  0j”.  As  a  final 
step,  we  make  sure  that  the  interpolated  longitude  in  the  double -primed  coordinate  system  0”  is 
between  -180°  and  +180°.  If  it  is  not,  then  360°  is  either  added  or  subtracted  to  the  value  of  0” 
until  it  is. 


Now  that  we  have  the  interpolated  point  (0,0”)  in  the  double -primed  coordinate  system 
(remember,  the  point  is  on  the  double-primed  equator),  we  reverse  the  two-step  process  to  obtain 
the  coordinates  (4>i, 9;)  in  the  original  coordinate  system.  We  shall  review  equations  158,  164, 
and  155.  Solving  equation  158  for  (j);,  we  see  that 

(j)i  =  arc  sin  (cos  0”  sin  4>max)  .  (173) 

The  range  of  (j);  is  between  -90°  and  +90°,  and  we  see  that  equation  173  delivers  a  unique  value 
for  4>i.  Now  that  we  know  the  value  of  cf);,  we  may  solve  for  9;.  Solving  equations  164  and  155 
for  0i,  we  find 


and 


0i  =  arc  sin  (sin  0;”  /  cos  cj)0  +  9o 


(174) 


0i  =  arc  cos  (cos  0i”  cos  cj)max  /  cos  cj); )  +  9o . 


(175) 
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Once  again,  we  find  that  9;  ranges  from  -180°  to  +  180°.  But  the  range  of  values  offered  by 
computer  languages  (29-32)  returns  values  between  -90°  and  +90°  for  arc  sin  and  0°  and  180° 
for  arc  cos.  Reusing  the  process  when  we  transformed  to  the  double-primed  coordinate  system, 
we  arrive  at  the  value  of  0;: 

if  sin  9”  /  cos  cj);  >  0,  then 
0i  =  arc  cos  (cos  9;”  cos  4>max  /  cos  ());)  +  0o , 
and  if  sin  9;”  /  cos  cj);  <  0,  then 

0;  =  -arc  cos  (cos  0”  cos  (t)max  / cos  cj); )  +  90 .  (176) 

Keep  in  mind  that  9;  must  be  between  -180°  and  +180°.  If  it  is  not,  then  add  or  subtract  360°  as 
is  appropriate.  We  have  finished  describing  the  process  of  obtaining  interpolated  points  along  a 
great  circle  on  a  sphere. 

4.3.5  Summary  of  Inputting  Waypoint  Data  and  Obtaining  Interpolated  Points 

This  section  briefly  describes  the  process  of  inputting  waypoint  data  and  obtaining  the  location 
of  the  interpolated  points  along  a  great  circle. 

Step  1 :  Assign  the  waypoint  longitudes  0j,  latitudes  cj)j,  and  elevations  E,.  If  the  waypoints  are 
for  a  ground  platform,  the  elevations  are  read  from  a  database. 

Step  2:  Calculate  the  great  circle  equation  that  connects  consecutive  waypoints.  Most  of  the 
time,  it  will  be  possible  to  use  equation  10  to  get  90  and  equation  1 1  to  get  (t)max,  which  are  the 
descriptive  parameters  of  equation  6,  the  great  circle  equation. 

Step  3:  Calculate  the  direct  distance  A  with  equation  79,  the  angular  distance  A0  with  equation 
80,  and  the  constant  k  with  equation  95  for  consecutive  waypoints.  This  makes  it  possible  to 
calculate  the  distance  Dj  between  consecutive  waypoints  with  equation  98. 

Step  4:  Obtain  the  cruising  velocity  vcru  and  acceleration  constant  ao  for  the  platfonn  from  a 
database.  Then,  assign  either  the  arrival  time  taj+i  to  each  waypoint  to  calculate  the  waypoint 
velocity  Vj+i  (which  can  usually  be  done  with  equation  76)  or  assign  the  waypoint  velocity  and 
calculate  the  arrival  time  taj+i  (which  can  usually  be  done  with  equation  75). 

Step  5:  Calculate  the  inflection  times  fi  and  t2,  which  can  usually  be  found  using  equations  72 
and  73. 

Step  6:  Select  the  loiter  time  (tdj  -  taj)  for  each  waypoint.  Notice  that  this  can  only  be  done  if 
the  waypoint  velocity  Vj  is  zero.  If  the  waypoint  velocity  Vj  is  not  zero,  the  loiter  time  defaults  to 
zero. 
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Note  that  steps  1-6  are  executed  every  time  a  waypoint  is  added  to  a  platform  path.  After  all  the 
waypoints  for  a  platfonn  have  been  entered,  proceed  with  step  7-10  and  calculate  the 
interpolated  latitudes  and  longitudes  (cj)i,0i). 

Step  7:  Select  the  interpolation  times,  (each  one  delineated  t,)  for  the  platform.  Use  these  times 
to  calculate  the  interpolated  distances  D;  with  equations  133-135.  Obtain  the  interpolated 
elevation  E;  with  equation  137  and  the  interpolated  angle  A0,  with  equation  139. 

Step  8:  Transform  the  waypoints  with  departure  times  tdj  and  arrival  times  taj+i,  which  bound 
the  interpolated  time  t,  to  a  new,  double-primed  coordinate  system  with  longitudes  0/’  and  0j+i” 
using  equations  165  and  166.  The  points  will  be  on  the  double-primed  equator. 

Step  9:  Calculate  the  interpolated  longitude  0;”  using  equation  172. 

Step  10:  Transform  the  interpolated  longitude  0;”  (it  too  is  on  the  double-primed  equator)  back 
to  the  original  coordinate  system  with  the  longitude  and  latitude  (4>i,0i)  using  equations  173  and 
176. 


5.  Conclusions 


The  algorithms  and  formulae  describing  the  motion  of  communications  platforms  needed  for  the 
NCAM  deployment  module  have  been  successfully  derived.  For  purposes  of  verification,  the 
details  of  the  algorithm  and  formulae  development  have  been  documented  and  are  shown  to  be 
consistent  with  the  assumptions  implicit  in  the  model’s  description  of  platform  motion.  The 
assumptions  were  that  the  platform  travels  along  great  circle  paths  between  waypoints,  the 
platforms  accelerate  and  decelerate  at  a  constant  rate  between  waypoints,  the  platforms  have  a 
maximum  cruising  speed  consistent  with  the  platform  and  terrain  type  (HMMWV  vs.  tank  vs. 
UAV  vs.  dismounted  Soldier  on  paved  road  vs.  rough  mountainous  terrain  vs.  flat  desert  terrain, 
etc.),  and  the  platforms  change  elevation  at  a  constant  rate. 

This  approach  allowed  the  user  to  outline  each  platfonn’s  path  with  waypoints  consistent  with 
the  user’s  specification  of  arrival  and  departure  times  at  each  waypoint.  Using  the  equations 
derived  provided  the  user  with  upper  and  lower  limits  for  inputted  data.  This  gave  the  user 
maximum  freedom  to  select  input  and  ensure  that  the  inputs  were  restricted  so  as  to  be  self- 
consistent.  This  made  it  possible  to  calculate  the  position  of  platforms  between  waypoints  at 
predetermined  snapshot  times.  Hence,  it  will  become  possible  to  model  link  viability  between 
moving  platforms  once  the  other  modules  are  completed.  Therefore,  the  equations  and 
algorithms  presented  here  are  mandatory  to  allowing  NCAM  to  model  communications  links 
between  moving  nodes. 
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Making  the  motion  model  more  complex  would  have  made  it  possible  to  model  platform  motion 
more  realistically.  However,  the  derivations  of  the  equations  and  algorithms  would  have  to  be 
far  more  involved  and  complex  than  the  ones  shown  here.  Furthermore,  it  would  have  demanded 
that  more  information  be  input  from  the  user.  Making  the  motion  simpler  would  have 
diminished  the  realism  of  the  platform  motion  and  diminished  the  control  the  user  has  over 
platform  motion.  The  approach  used  here  (specifying  waypoint  location  and  the 
arrival/departure  velocity  or  time)  provided  the  best  compromise  between  user  control,  equation 
and  algorithm  simplicity,  and  realistic  platform  motion. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


CCI 

connectivity  confidence  interval 

DTED 

Digital  Terrain  Elevation  Data 

FORTRAN 

Fonnula  Translator 

HMMWV 

high-mobility  multipurpose  wheeled  vehicle 

NCAM 

Network  Connectivity  Analysis  Model 

S/N 

signal  to  noise 

TIREM 

Terrain  Integrated  Rough  Earth  Model 

UAV 

unmanned  aerial  vehicle 

WGS84 

World  Geodetic  System  for  1984 
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